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ABSTRACT 


We  wish  to  solve  fluid  flow  problems  in  only  a  portion  of  a  large  or  infinite 
domain.  By  restricting  our  area  of  interest,  we  effectively  create  a  boundary  where 
none  exists  physically,  dividing  our  computational  domain  from  the  rest  of  the  physical 
domain.  The  challenge  we  must  overcome,  then,  is  defining  this  boundary  in  such  a 
way  that  it  behaves  computationally  as  if  there  were  no  physical  boundary.  Such  a 
boundary  definition  is  often  called  a  non-reflecting  boundary,  as  its  primary  function 
is  to  permit  wave  phenomena  to  pass  through  the  open  boundary  without  reflection. 
In  this  dissertation  we  develop  several  non-reflecting  boundary  conditions  for  the 
linearized  Euler  equations  of  inviscid  gas  dynamics.  These  boundary  conditions  are 
derived  from  the  Higdon,  Givoli-Neta,  and  Hagstrom-Warburton  boundary  schemes 
for  scalar  equations,  and  they  are  adapted  here  for  a  system  of  first-order  partial 
differential  equations.  Using  finite  difference  methods,  we  apply  the  various  boundary 
schemes  to  the  gas  dynamic  equations  in  two  dimensions,  in  an  open  domain  with 
and  without  the  influence  of  gravity  or  Coriolis  forces.  These  new  methods  provide 
significantly  greater  accuracy  than  the  classic  Sommerfeld  radiation  condition  with 
only  a  modest  increase  to  the  computation  time. 
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I.  INTRODUCTION 

Many  problems  in  compntational  flnid  dynamics  occnr  within  a  limited  portion 
of  a  very  large  or  infinite  domain.  Difficnlties  immediately  arise  when  one  attempts 
to  dehne  the  bonndary  condition  for  snch  a  system.  Snch  bonndary  conditions  are 
necessary  for  the  problem  to  be  well-posed,  bnt  the  physical  system  nnder  consider¬ 
ation  has  no  bonndary  to  model.  One  needs  to  dehne  an  artihcial  bonndary  whose 
behavior  models  the  open  edge  of  the  physical  system.  Snch  a  bonndary  dehnition  is 
often  called  a  non-reflecting  bonndary  condition  (NRBC),  as  its  primary  fnnction  is 
to  permit  wave  phenomena  to  pass  throngh  the  open  bonndary  withont  rehection. 

NRBC  development  has  been  an  ongoing  research  area  since  the  1960s.  As  with 
any  compntational  endeavor,  there  are  trade-offs  involved.  An  ideal  NRBC  wonld  be 
fast,  accnrate,  stable,  and  easy  to  implement:  fast,  meaning  that  the  compntation 
of  bonndary  valnes  is  small  or  negligible  relative  to  the  interior  domain;  accurate, 
meaning  that  there  is  little  to  no  spnrions  rehection  indnced  by  the  bonndary  con¬ 
dition;  stable,  meaning  that  the  bonndary  compntations  do  not  canse  the  solntion 
to  degrade  catastrophically  over  time;  and  easy  to  implement,  meaning  that  the  nser 
can  incorporate  the  NRBC  compntations  into  an  operational  model  with  minimal 
modihcation  to  existing  code.  Realistically,  one  mnst  settle  for  two  or  three  of  these 
attribntes,  at  best.  Conseqnently,  the  search  for  better  NRBC  methods  continnes.  In 
addition,  researchers  continne  to  apply  existing  NRBC  methods  to  new  domains  and 
wave  propagation  eqnations. 

In  this  dissertation  we  develop  several  NRBCs  for  the  linearized  Enler  eqna¬ 
tions  of  inviscid  gas  dynamics.  These  bonndary  conditions  are  derived  from  the 
Higdon,  Givoli-Neta,  and  Hagstrom-Warbnrton  bonndary  schemes  for  scalar  eqna¬ 
tions,  adapted  here  for  a  system  of  hrst-order  partial  diherential  eqnations.  Using 
hnite  diherence  methods,  we  apply  the  varions  bonndary  schemes  to  the  gas  dynamic 
eqnations  in  two  dimensions,  in  an  open  domain  with  and  withont  the  inflnence  of 
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gravity  or  Coriolis  forces.  These  new  methods  provide  signihcantly  greater  accuracy 
than  the  classic  Sommerfeld  radiation  condition  with  only  a  modest  increase  to  the 
computational  cost. 

Our  motivation  for  developing  these  NRBCs  is  to  support  the  work  of  Giraldo 
and  Restelli  in  their  efforts  to  develop  the  next  generation  of  mesoscale  atmospheric 
modeling  tools  [32].  In  their  current  form,  the  models  rely  on  large  absorbing  layers 
surrounding  the  computational  domain.  In  order  to  be  effective,  these  “sponge  lay¬ 
ers”  can  be  as  thick  as  the  original  domain  [31],  resulting  in  a  total  domain  which 
is  nearly  quadrupled  in  size.  If  we  can  replace  these  large  sponge  layers  with  accu¬ 
rate  NRBCs,  then  the  modeling  tools  will  requires  less  memory  and  computational 
overhead,  signihcantly  increasing  their  efficiency. 

The  rest  of  the  dissertation  is  outlined  as  follows.  We  begin  in  Chapter  II  by 
deriving  the  equations  under  consideration,  the  linearized  Euler  equations  of  inviscid 
gas  dynamics.  In  Chapter  III  we  provide  an  overview  of  existing  NRBCs  and  provide 
specihc  details  about  the  Higdon,  Givoli-Neta,  and  Hagstrom-Warburton  schemes  for 
scalar  equations.  We  then  extend  these  boundary  conditions  to  the  hrst-order  2-D 
linearized  Euler  system  in  Chapters  IV  (Higdon),  V  (Givoli-Neta),  and  VI  (Hagstrom- 
Warburton).  In  all  three  cases,  we  consider  the  NRBC  implementation  in  a  semi- 
inhnite  or  inhnite  channel  and  in  an  open  domain,  under  basic  conditions  as  well 
as  under  the  inhuence  of  Coriolis  forces  or  gravity.  Numerical  examples  using  hnite 
differences  are  provided  throughout.  We  discuss  the  issue  of  long-time  stability  in 
Chapter  VH.  We  offer  a  qualitative  comparison  of  the  three  NRBC  techniques  in 
Chapter  VHI.  We  close  in  Chapter  IX  with  a  summary  of  our  results  and  a  list  of 
areas  for  further  research. 
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II. 


MODELING  INVISCID  FLUID  FLOW 


In  this  chapter  we  explore  the  equations  governing  the  motion  of  a  body  of 
fluid.  These  principles  describe  the  flow  of  water  in  a  channel  or  in  the  ocean,  air 
movement  over  mountains,  airflow  and  drag  in  aircraft  design,  and  even  the  heat 
generated  by  a  spacecraft  re-entering  the  atmosphere.  Although  many  physical  phe¬ 
nomena  depend  on  the  viscosity  of  the  fluid,  certain  large-scale  flows  of  low- viscosity 
fluids  (e.g.,  air)  can  be  reasonably  approximated  by  assuming  the  viscosity  is  negli¬ 
gible. 

By  neglecting  viscous  forces,  our  fluid  flow  equations  can  be  derived  based 
simply  on  physical  conservation  laws  governing  mass,  momentum,  and  energy.  The 
following  section  considers  each  conservation  law  in  turn  and  derives  the  relevant 
governing  equations  therefrom. 

We  derive  the  Euler  equations  based  hrst  on  internal  factors.  Then  we  consider 
the  inhomogeneous  factors  which  affect  these  equations  in  the  context  of  atmospheric 
modeling. 


A.  DERIVATION  OF  EQUATIONS 

1.  Conservation  of  Mass 


Mass  is  conserved  in  a  closed  system.  If  dm  denotes  an  inhnitesimal  portion 
of  the  mass,  then  /  dm  denotes  the  total  mass  of  the  body,  and  conservation  of  mass 
requires 


d 

dt 


dm  =  0 


(11.1) 


Since  the  mass  is  continuous,  we  can  bring  the  derivative  inside  the  integral;  thus 


(11.2) 
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Furthermore,  since  this  statement  must  be  true  for  any  piece  of  the  body’s  mass,  the 
integrand  must  be  identically  zero,  i.e.. 


or 


d  ,  ,  X 
—  (dm)  =  0 
dV  ^ 


—  {pdxdydz)  =  0 

CLL 


Applying  the  product  rule  gives 


f  1  1  1  \dp  /  7  7  X  ,d^  /  7  7  \  Ay  /  7  7  \  A^ 

(dxdydz)—  +  (pdydz)d—  +  (pdxdz)d—  +  (pdxdy)d—  =  0 

Dehne  the  velocity  vector  u  =  (u,  n,  where  u  =  dx/dt,  v  =  dy/dt,  and  w 
and  divide  all  four  terms  by  dxdydz  to  get 

dp  du  dv  dw 

Applying  the  chain  rule  to  each  derivative  gives 

7^  +  u^  +  v^  +  wlf  +  p(P'A  +  P^  +  7r^)  ' 

at  ox  oy  oz  •  \ox  ax  oy  dx  oz  dx } 

-Ln  ( \  ^  ( I  dtv^  i  dw  dz\ 

\dx  dy  dy  dy  '  dz  dy )  '  r  \  dx  dz  '  dy  dz  '  dz  dz )  ) 

However,  since  x,  y,  and  ^  are  independent,  this  equation  reduces  to 
dtp  +  udxp  +  vdyp  +  wdzp  +  p  {dxU  +  dyV  +  d^w)  =  0, 


y  =0 


i.e., 


or,  in  vector  notation. 


OtP dy{fm}  +  ^  0 


ftp  +  V  ■  (pm)  =  0  , 


(11.3) 

(11.4) 

(11.5) 
dz/dt, 

(11.6) 

(11.7) 

(11.8) 

(11,9) 

(11,10) 


where  from  here  on,  we  use  the  following  shorthand  for  partial  derivatives: 

f,  _8_  a 

*  dt  '  dxdy 

This  equation  assumes  no  external  sources  or  sinks  adding  or  removing  mass  from 
the  system.  We  note  in  passing  that  we  have  shown,  from  first  principles,  a  special 
case  of  Reynold’s  transport  theorem  [3]  applied  to  the  density  of  a  fluid. 
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2.  Conservation  of  Momentum 

We  now  consider  the  momentnm  of  the  flnid  body.  To  this  end,  we  hrst 
determine  what  forces  act  npon  the  flnid.  For  the  pnrpose  of  this  dissertation,  the 
flnid  body  is  assnmed  to  be  a  portion  of  the  Earth’s  atmosphere. 

a.  Forces  Acting  upon  a  Fluid  Body 

The  internal  forces  exerted  by  the  body  on  itself  consist  of  pressnre 
forces  and  viscons  forces  (which  are  neglected  since  we  are  assnming  an  invsicid 
flnid).  The  external  inhomogeneons  forces  we  consider  are  the  effects  of  gravity  and 
the  Earth’s  rotation. 

(1)  Force  Dne  to  Internal  Pressnre.  The  pressnre  force 
resnlts  from  pressnre  differences  within  the  body  (see  Fig.  1)  and  acts  to  retard  the 
motion  of  the  flnid. 


Fignre  1.  Pressure  Differences  in  a  Small  Volume  [From  [113],  Fig.  3,  p.  15] 

In  a  small  piece  of  fluid  volume  dV  =  dxdydz,  the  change  in 
pressure  in  the  x  direction  is  (to  hrst  order  accuracy) 

dxP  dx 
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The  force  exerted  by  this  pressure  difference  is  the  product  of  the  pressure  difference 
and  the  surface  area  on  which  the  pressure  is  exerted.  Hence,  the  pressure  force  in 
the  X  direction  is  given  by 


dFpressurex  dx^  dydz  ^xP  dV 

Similarly,  the  forces  exerted  by  pressure  differences  in  the  y  and  ^  directions  are  given 
by 

dFpressurey  =  “  {dyP  dy)  dxdz  =  -dyP  dV 

dFpxessurcz  i^^zP  dz^  dxdy  dzP  dV  ^11.12^ 


Thus,  the  total  force  due  to  internal  pressure  differences  is  the  sum  of  the  three 
components. 


dFpressure  =  “  (dxPl  +  OyPj  +  O^pk^  dV  =  -{Vp)dV, 


and  the  overall  force  on  the  entire  body  is 


pressure 


J  dFpressure  J  {^p^dV 


(11.13) 


(11.14) 


(2)  Force  Due  to  Gravity.  According  to  Vallado  [112],  the 
gravitational  acceleration  comes  from  Newton’s  Law  of  Gravitation: 


a  = 


/i 


r 

If  I 


(11.16) 


where  p  is  the  Earth’s  gravitational  parameter  (398,600.4418  km^/s^),  and  ||f||  is  the 
radius  from  the  center  of  the  Earth  (6378.137  km  at  sea  level).  At  sea  level,  this  gives 
a  gravitational  acceleration  of  approximately  9.798  m/s^.  Gomparing  the  magnitudes 
of  the  acceleration  at  sea  level  to  that  at  an  altitude  of  20  km,  we  find 


Aa  =  /i 


1 


1 


(r  +  20km)  2 


m 

-0.06116  — 


(11.16) 


Given  such  a  small  relative  difference,  we  can  treat  the  gravitational  acceleration  as 
constant,  so  the  gravitational  force  is  simply 


gravity 


-gdmk 


(11.17) 
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where 


g  ^  9.798—  (11.18) 

s 

(3)  Force  Due  to  Earth’s  Rotation.  In  Fig.  2,  D  repre¬ 
sents  the  angular  velocity  of  the  sphere,  which  for  Earth  has  a  value  of  27r  radians  per 
sidereal  day  (23  hours,  56  minutes,  4.090524  seconds  [112])  or  7.292116  x  10“®  rad/s. 
We  must  consider  the  effects  caused  by  the  rotating  reference  frame.  Here  we  follow 


Figure  2.  Rotating  Sphere  [From  [113],  Fig.  2,  p.  11] 

the  derivation  laid  out  by  Holton  [70]  and  Pedlosky  [94].  Let  r  denote  the  position 
vector  of  an  element  in  our  rotating  reference  frame  (see  Fig.  3,  where  'ijj  is  the  angle 
between  Fand  the  angular  rotation  vector  D).  In  a  small  time  interval  At,  the  vector 
rotates  by  an  angle  A6  =  |D|At.  If  the  time  interval  is  small  enough,  the  change  in 
the  vector  is  given  by 

Af=  n\f\A6  sin 'ijj , 

where  n  is  the  unit  vector  in  the  direction  of  the  cross  product  D  x  r ,  i.e., 

^  Q  X  f 

n  =  - 

|D  X  r] 


(11.19) 

(11.20) 
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n 


Figure  3.  Vector  in  Rotating  Reference  Frame  [After  [94],  Fig.  1.5.2,  p.  15] 


Thus, 


or  in  the  limit  as  At  — )■  0, 


Ar 


Q  X  f 
X  r] 


rl-— sin-iw, 
'  'At 


df 

dt 


Q  X  f 
X  r] 


|r1 


dO 


(11.21) 


(11.22) 


and  since  ^  =  |fl|,  we  get 


dr  ^  ^ 

—  =  Vt  X  r 
dt 


(11.23) 


If  we  have  f  =  xt  +  yj  +  zk  in  the  east-north-up  coordinate  frame  shown  in  Fig.  2, 
then  we  convert  from  rotating  to  inertial  reference  frames  by 


^  dr  ^  ^ 

ui  =  u-\-  —  =  u  +  \lxr^ 
dt 


(11.24) 


where  u  is  the  change  in  position  within  the  rotating  reference  frame,  and  uj  is  the 
same  motion  relative  to  the  inertial  reference  frame.  Differentiation  of  the  above. 
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substituting  uj  for  r  and  noting  that  is  constant,  yields 


dui\  / dui\  -  ^ 

(j  ‘7  /  — ^  ^  ^ 

=  —  +  2nxu  +  nx{nxf^ 

CLL 


(11.25) 


Based  on  Fig.  2,  our  angular  momentum  vector  is 


12  =  1 12 1  cos  (pj  +  1 12 1  sin  (f)k 


(11.26) 


Using  this  dehnition  and  the  dehnition  of  the  cross  product  twice,  we  have 

/  rr  \ 

n  X  (12  X  r)  =  |12|^  ^  sin  0  cos  0  —  1/ sin^  0  ,  (11.27) 

\  — ^  cos^  (f)  +  ysincf)  cos  (f)  / 

but,  since  |12p  =  5.317496  X  10  we  can  ignore  the  last  term  of  (11.25)  to  get 


/  dui\  du 

-r-  =  —  +  212  X  n 

\  dt  J  j  dt 

Applying  this  formula  to  the  momentum  vector  gives 


(11.28) 


=  ^{udm)  +  212  X  (udm)  (11.29) 

CLL 

Hence,  when  we  consider  the  change  in  momentum,  we  must  add  to  it  this  rotational 
effect.  Now  this  rotational  effect  can  be  simplihed.  The  atmosphere  is  thin  compared 
to  the  radius  of  the  Earth.  Furthermore,  atmospheric  flows  tend  to  be  parallel  to  the 
ground  (i.e.,  the  velocity  in  the  ^  direction  is  small).  Taking  the  cross  product  of  this 
vector  and  our  velocity  vector  gives 


2(u,drn) 


212  X  -u  =  2|12|((wcos0  —  usin0)£  +  usm(f)j  —  u  cos  (j)k) 


(11.30) 


Again,  assuming  a  thin  atmosphere  and  thus  neglecting  terms  in  the  ^  direction,  this 
simplihes  to 

2Cl  X  u  ^  2|12|  sm(p{—vi  +  uj)  =  2|12|  sm(j{k  x  u)  (11.31) 
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If  we  define  /  =  2|fl|  sin0,  we  can  express  this  force  as 


d 


dt 


{udm)\  +  f{Iv  X  (udm)) 


7 


(11.32) 


R 


Since  we  are  observing  the  flow  from  the  rotating  reference  frame,  this  additional  force 
becomes  an  inhomogeneons  term  on  the  right-hand-side  of  the  momentnm  eqnation: 


Fcorioiis  =  J  f{ux  k)dm 

b.  Summary  of  Forces 

Thns,  the  eqnation  for  conservation  of  momentnm  is 


(11.33) 


7 

dt 


J  udm  =  —  J  WpdV  -|-  j  f{u  x  k)dm  —  J 


gkdm 


Here,  as  before,  the  integrand  mnst  be  identically  zero. 


^iudm)  -f  ^^dm  =  fiu  x  k)dm  —  gkdm, 
dt  p 


(11.34) 


(11.35) 


or,  by  component. 


^(udm)  -f  =  fvdm 

dt  p 

d  din 

—  {vdm)  H - dyp  =  —fudm  (11.36) 

d  dm  ^ 

—  iwdm)  - OzP  =  —gdm 

dt  p 

Using  a  derivation  similar  to  that  in  Section  1,  these  eqnations  can  be  rewritten  in 
vector  notation: 


dt{pu)  +  V  •  {pu  <^u)  +  Vp  =  f{pu  X  k)  —  gpk, 


(11.37) 


where  0  denotes  the  tensor  prodnct.  However,  we  can  simplify  the  eqnations  to 


dtu  +  udxU  +  vdyU  +  wdzU  H — d^p  =  fv 

P 

dtv  +  udxV  -|-  vdyV  +  wdzV  -I — dyP  =  —fu 

dtw  +  udxW  -|-  vdyW  -f-  wdzW  H — dzP  =  —g, 

P 


(11.38) 


or,  in  vector  notation. 


du  1  /V  /V 

—  -1-  (m  •  V)m  -f  -Vp  =  fiu  X  k)  —  gk 

dt  p 


(11.39) 
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3.  Conservation  of  Energy 

The  intrinsic  energy  in  each  piece  of  mass  dm  consists  of  its  kinetic  energy.  For 
atmospheric  modeling,  there  are  also  inhomogeneons  components  for  thermal  energy 
and  potential  energy.  These  qnantities  are  dehned  by 

dHE  =  CyTdm 

IliTIP 

dKE  =  (11.40) 

Li 

dPE  =  gzdm 

In  these  eqnations,  Cy  is  the  specihc  heat  at  constant  volnme,  and  T  is  the  temperatnre. 
Hence  the  total  energy  in  the  body  is 

J  {^yT  +  dm 

To  change  the  amonnt  of  energy,  apply  force  to  the  body  over  a  distance  or  add/remove 
heat.  Adding/removing  heat  is  a  complicated  process  involving  solar  radiation,  ther¬ 
mal  radiation  from  the  Earth,  heat  dissipation  into  space,  reflectivity  of  the  Earth’s 
snrface,  and  other  factors.  Dne  to  the  complexity,  factors  which  contribnte  to  a  change 
in  temperatnre  will  not  be  considered  here.  Rather,  we  will  simply  consider  the  time 
derivative  of  onr  heat  energy  term  on  its  own  withont  determining  the  precise  factors 
which  determine  it.  Later  in  this  section,  we  will  hide  the  temperatnre  component 
entirely,  combining  thermal  and  kinetic  energy  into  a  “total  energy”  term  e. 

For  force  application,  the  momentnm  eqnations  assnme  a  balance  between 
the  external  forces  and  the  internal  pressnre.  However,  if  the  two  qnantities  are  not 
matched,  then  the  volnme  will  expand/contract  so  that  they  are  matched.  If  there 
is  an  imbalance  between  the  forces  on  opposing  sides  of  the  body  (hence  a  pressnre 
difference),  then  the  body  will  accelerate  in  the  direction  of  the  net  force.  (Imagine 
pnlling  a  spring  across  a  table.  The  spring  will  stretch  according  to  the  pnlling  force, 
and  the  stretched  spring  will  move  in  the  pnlling  direction.) 

We  can  conceive  of  the  total  force  as  having  these  two  components,  compression 
and  acceleration.  For  compression  at  constant  temperatnre,  as  the  volnme  decreases. 
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the  internal  pressnre  increases,  and  vice  versa.  Hence, 


(11.41) 


However,  if  the  temperatnre  is  not  constant,  then  the  pressnre  increases  proportionally 
to  the  thermal  energy  at  the  constant  volnme,  adding  the  following  component  to  onr 


eqnations: 


^(c^Tdm)  =  ^dV 
dV  ^  dt 


(11.42) 


The  acceleration  component  is  determined  by  Newton’s  Second  Law.  The  acceleration 
occnrs  as  a  resnlt  of  an  imbalance  between  the  force  acting  on  one  side  of  the  body 
and  the  force  acting  on  the  opposite  side,  which  resnlts  in  the  pressnre  differences 
dehned  in  the  conservation  of  momentnm  section.  Using  F  =  ma  and  onr  dehnition 


of  dFpressure  WC  have 


-{'Vp)dV  Vp 

dm  p 


(11.43) 


This  force  changes  the  kinetic  energy  of  the  body.  The  time  derivative  of  kinetic 


energy  is 


—  {KE)  =  dm  iu  ■  =  dm{u  ■  a) 


Using  onr  formnla  for  the  acceleration,  we  have 


j^{KE)  =  -^{u  ■Vp)  =  -{u-  Vp)dV 


(11.44) 


(11.45) 


In  the  context  of  atmospheric  modeling,  this  force  also  resnlts  in  an  increase  in  the 


body’s  potential  energy  as  it  changes  position  vertically,  so  that  we  have 

^^{KE)  +  j^{PE)  =  -{u-Vp)dV 


(11.46) 


Therefore,  changing  the  energy  in  the  body  resnlts  in  increased  kinetic  energy,  body 
compression,  and  pressnre  increase  dne  to  temperatnre 


_  /  ^  +  sd  dmU  / -ipiV)  =  -  /(«.  Vp)dV  +  /  liV  (11.47) 
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Combining  the  integrals  and  bringing  the  derivative  inside  gives 


/ 1  + 

where  we  use  the  product  rule  to  combine 
integral  must  be  true  for  any  body,  so  the 

d((c„r+l^+9.)dm) 


{u  ■Vp)dV  +  p-^{dV)  =  0,  (11.48) 

^{pdV)  —  ^dV  into  p-^{dV).  Again,  this 
integrand  must  be  identically  zero 

+  {u-Vp)dV +  p^{dV)  =  Q  (11.49) 

CLL 


II  “'ll  2 

Let  e  =  CvT  +  denote  the  total  internal  energy  of  the  system  (i.e.,  not  including 
potential  energy).  Then  we  have 


A 

dt 


((e  +  gz)pdV)  +  {ud^p  +  vdyp  +  wdzp)dV  +p^{dV) 

CLL 


0 


(11.50) 


It  is  easy  to  show  that  this  equation  is  equivalent  to 


dt{pe)  +  dx{{pe  +  p)u)  +  dy{{pe  +  p)v)  +  dz{{pe  +  p)w)  =  -gpw,  (11.51) 
or,  in  vector  notation, 

dt{pe)  +  V  ■  ((pe  +  p)u)  =  -(gk)  ■  (pu)  (11.52) 

This  equation  can  be  placed  in  a  simpler  form,  if  we  assume  our  fluid  is  an  ideal 
gas.  Using  the  ideal  gas  law  p  =  pRT  [24,  58],  where  R  =  Cp  —  c^,  simple  algebraic 
manipulation  using  the  previous  equations  can  reduce  this  energy  equation  to 


dtp  +  udxP  +  vdyp  +  wdzP  +  7P  {dxU  +  dyV  +  dzw)  =  0  ,  (11.53) 

where  7  =  — .  In  vector  notation,  this  equation  can  be  written  as 

Cy 

dtp  +  {u  ■  V)p  +  7p(V  •  -u)  =  0  (11.54) 

While  this  equation  is  simple,  it  lacks  an  explicit  energy  term  and  thus  fails  to  convey 
the  energy  conservation  principle  from  which  it  was  derived. 
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B.  SUMMARY  OF  NON-LINEAR  EULER  EQUATIONS 

Using  only  the  physical  conservation  principles  for  mass,  momentnm,  and  en¬ 
ergy,  we  can  derive  Enler’s  eqnations  for  inviscid  flnid  motion.  We  have  a  system  of 
hve  eqnations  for  six  variables: 

dtp  +  da;{pu)  -1-  dy{pv)  -f  d^{pw)  =  0 

dt{pu)  +  d^{pu^)  -f  dy{puv)  +  d^{puw)  +  d^p  =  fpv 
dt{pv)  +  d^{puv)  +  dy{pv‘^)  +  d^{pvw)  +  dyp  =  -fpu  (11.55) 

dt{pw)  +  dx{puw)  +  dy{pvw)  +  dz{pw^)  d^p  =  -gp 
dt{pe)  +  da;{{pe  +  p)u)  +  dy{{pe  +  p)v)  +  dz{{pe  +  p)w)  =  -gpw 

In  vector  notation,  the  eqnations  can  be  written  as 

dtp  +  V  ■  {pu)  =  0 

dt{pu)  +  V  ■  {pu  ®u)  +  Vp  =  f{pu  X  k)  —  gpk  (11.56) 

dt{pe)  +  V  ■{{pe  +  p)u)  =  -{gk)  ■  {pu) 

A  state  eqnation  of  some  kind  is  reqnired  to  close  the  system.  For  an  ideal  gas,  we 
can  nse  p  =  pRT  to  simplify  the  eqnations  to 

dtp  +  d^{pu)  +  dy{pv)  +  d^{pw)  =  0 

dtu  -f  udxU  +  vdyU  +  wdzU  -\ — dxP  =  fv 

P 

dtv  +  udxV  +  vdyV  +  wdzV  H — dyP  =  —fu  (11.57) 

dtw  +  udxW  +  vdyW  +  wdzW  H — dzP  =  —g 

P 

dtp  +  udxP  +  vdyp  +  wdzP  +  7P  {dxU  +  dyV  +  dzw)  =  0 

(see  Appendix  A  for  details).  For  all  three  formnlations,  the  terms  on  the  right-hand- 
side  denote  forces  specihc  to  atmospheric  modeling. 

In  the  next  section,  we  will  derive  the  linearized  form  of  (11.57),  which  will 
form  the  basis  for  developing  the  hnite-difference  implementation  of  the  NRBCs  in 
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Chapters  IV-VI.  We  will  begin  with  a  simplified  prototype  implementation  of  the 
antomated  Higdon  NRBCs  with  no  advection  or  forcing  terms,  and  we  will  bnild 
the  implementation  to  inclnde  Coriolis,  gravity,  and  non-zero  mean  flow.  We  will 
also  develop  anxiliary  variable  forms  which  eliminate  the  high-order  derivative  terms, 
nsing  both  the  Givoli-Neta  and  the  Hagstrom-Warbnrton  anxiliary  variable  NRBC 
formnlations. 

C.  LINEARIZED  EULER  EQUATIONS 

Having  derived  the  non-linear  eqnations,  we  now  seek  to  simplify  them  into  a 
linear  form.  We  assnme  that  each  state  variable  consists  of  a  time-independent  mean 
valne  and  a  small  pertnrbation  from  that  mean.  Thns, 

ip  =  (p  +  ip*,  (11.58) 

where  the  overbar  denotes  the  reference  valne,  and  the  asterisk  denotes  the  0(S) 
pertnrbation  variable.  Before  we  can  derive  this  linearized  form,  we  mnst  first  define 
onr  reference  variables.  We  will  perform  the  linearization  on  (11.57). 

1.  Defining  the  Reference  Variables 

The  reference  valnes  are  time-independent  by  definition,  bnt  they  may  not 
necessarily  be  constant  in  space.  It  is  reasonable  to  believe  that  the  reference  valnes 
for  the  velocity  variables  u,  v,  and  w  will  be  constant;  however,  this  may  not  be  trne 
for  the  density  p  and  pressnre  p.  In  the  presence  of  gravity,  a  volnme  of  compressible 
finid  will  be  compressed  by  the  weight  of  the  finid  above  it,  increasing  the  density  and 
pressnre.  Therefore,  it  is  reasonable  to  expect  that  onr  reference  states  for  density 
and  pressnre  will  vary  with 

Let  ns  consider  a  constant  domain  governed  by  (11.57).  With  everything  at 
rest,  we  set  all  time  derivative  terms  and  all  velocity  terms  to  zero.  This  leaves  ns 
with 


dxP  =  0 
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dyP  =  0 


(11.59) 


dzP  =  -gp 

This  confirms  our  expectation  that  our  reference  values  for  p  and  p  are  dependent 
on  So  for  our  reference  values,  we  let  mq,  vq,  and  wq  define  our  constant  mean 
velocities;  p  denotes  our  ^-dependent  density  reference  state,  with  po  the  density  at 
^  =  0;  and  p  denotes  our  ^-dependent  pressure  reference  state,  with  po  the  pressure 
at  ^  =  0. 

2.  Linearizing  the  Equations 

Begin  with  the  hrst  equation  of  (11.57): 

dtp  +  dx{pu)  +  dy{pv)  +  d:,{pw)  =  0.  (11.60) 

Expand  the  derivatives  via  the  product  rule,  and  substitute  the  reference/perturbation 
variable  expansion  in  (11.58)  to  get 

(p  +  P* )  +  (p  +  P* )  dx {uq  +  )  +  {uq  +  u*)  dx{p  +  p*) 

+  {p  +  p*)  dy{vo  +  V*)  +  {vo  +  V*)  dy{p  +  p*)  >=0.  (11.61) 

+  (p  +  p*)  d^{wo  +  w*)  +  {wo  +  w*)  d^{p  +  p*) 

/ 

Recalling  which  reference  variables  are  independent  in  space  and  time,  and  eliminating 
all  terms  of  0(5^),  we  get 

dtp*  +  P  {dxU*  +  dyV*  +  dz.w*)  +  u^dxP*  +  v^dyp*  +  w^d^p*  =  -d^p  {wq  +  w*) .  (11.62) 

We  leave  the  equation  in  this  form,  rather  than  reverting  it  to  the  whole  state  variable. 
By  considering  only  the  perturbation  variable,  we  eliminate  the  possibility  of  the 
reference  value  overwhelming  the  perturbation,  introducing  unnecessary  round-off 
errors  into  the  hnite- precision  calculations. 

When  we  apply  this  same  process  to  the  velocity  and  pressure  equations  of 
(11.57),  using  (11.59)  to  simplify  the  equation  for  w,  we  get 

dt<p  +  AdxP  +  Bdy(p  +  CdzP  =  D(p  +  E,  (11.63) 
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where 


A 


B 


C 


D 


E 


P*  U*  V*  w*  p*  ^ 

^  uq  p  0  0  0  ^ 

0  Mo  0  0  i 

0  0  Mo  0  0 

0  0  0  Mo  0 

^  0  7P  0  0  Mo  y 

^  Mo  0  p  0  0  ^ 

0  Mo  0  0  0 

0  0  Mo  0  4 

u  p 

0  0  0  Mo  0 

^  0  0  7P  0  Mo  y 

^  tMo  0  0  p  0  ^ 

0  tMo  0  0  0 

0  0  tMo  0  0 

0  0  0  tMo  4 

^  0  0  0  7P  tMo  y 

^  0  0  0  -d,p  0  ^ 

0  0/00 
0-/0  0  0 

-c/4  0  0  0  0 

^  0  0  0  0  y 

^  -dzpwo  ^ 

fvo 

-fuQ 

0 

^  gPWo  y 


(11.64) 
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For  the  ^{d^p,  dyp,  d^p}  terms  in  the  momentum  equations,  we  use  the  following 
linearization: 


dxP 

P 


P 


dzP 

P 


dx  {p  +  P* 
p  +  p* 
dxP* 


P  \1  +  ^ 


1--+  -  - 
P  \P 


dxP* 

P 

dxP* 

P 

dyP* 

p 


d,p*  +  d,p 

p  V  p  VP, 

dzP*  9P  L  P*\ 


P 


P 


P 


dzP*  I,  P*\ 


where  we  use  (11.59)  in  the  next  to  last  line. 

Note  that  the  matrix  A  is  singular  if  mq  =  0  or  Mq  =  7p/p.  Similarly,  B  is 
singular  if  uq  =  0  or  =  7p/p,  and  C  is  singular  if  wq  =  0  or  =  7p/p. 

So  what  have  we  lost  by  this  linearization?  The  main  difference  between  non¬ 
linear  flow  and  linear  flow  is  the  non-linear  interaction  between  vortices  [69].  For 
example.  Fig.  4  shows  a  rising  thermal  bubble  using  non-linear  equations  (left)  and 
their  linearized  form  (right).  The  turbulence  is  clearly  absent  in  the  linearized  case. 
However,  wave  motion  is  not  noticeably  affected.  Fig.  5  shows  an  acoustic  wave  using 
the  non-linear  (left)  and  linearized  equations  (right).  The  differences  are  negligible. 
(These  hgures  were  generated  during  some  early  work,  applying  a  spectral  element 
implementation  of  the  non-linear  system  (11.55)  in  the  xz  plane  influenced  by  gravity, 
and  to  the  linearized  form  of  the  same  system.)  Since  our  non-reflecting  boundary 
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conditions  are  intended  for  wave  propagation,  we  do  not  need  to  keep  the  non-linear 
effects  of  vorticity. 


Fignre  4.  A  Rising  Thermal  Bnbble  Using  Non-linear  (left)  and  Linear  (right)  Eqna- 
tions 


Fignre  5.  An  Aconstic  Wave  Using  Non-linear  (left)  and  Linear  (right)  Eqnations 

Before  we  begin  developing  NRBCs  for  this  eqnation  system,  let  ns  spend  some 
time  discnssing  NRBCs  in  general,  with  particnlar  emphasis  on  the  scalar-eqnation 
implementations  of  the  NRBCs  we  are  here  adapting  for  a  linear  hrst-order  system. 
This  discnssion  will  be  the  snbject  of  the  next  chapter,  and  then  the  new  NRBC 
implementations  will  be  developed  in  Chapters  IV-VI. 
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III.  NON-REFLECTING  BOUNDARY 

CONDITIONS 

A.  OVERVIEW  AND  HISTORY 

We  wish  to  solve  fluid  flow  problems  in  only  a  portion  of  a  large  or  infinite 
domain.  By  restricting  our  area  of  interest,  we  effectively  create  a  boundary  B  where 
none  exists  physically,  dividing  our  computational  domain  from  the  rest  of  the 
domain  T).  Implicit  in  our  choice  of  B  is  the  assumption  that  everything  we  wish  to 
model  is  contained  inside  fl;  nothing  external  impinges  on  our  computational  domain. 
The  challenge  we  must  overcome,  then,  is  defining  B  in  such  a  way  that  it  behaves 
computationally  as  if  there  were  no  physical  boundary.  How,  then,  do  we  specify 
what  occurs  at  these  boundaries?  The  usual  answer  is  to  claim  that  waves  flow  out 
of  the  domain  at  the  boundary,  but  they  do  not  flow  into  the  domain.  If  defined 
correctly,  this  claim  will  result  in  waves  which  propagate  out  of  the  computational 
domain  without  any  reflection,  so  that  the  computational  boundary  is  transparent 
to  these  outgoing  waves,  mimicking  the  real-world  behavior  where  no  such  physical 
boundary  exists. 

In  general,  there  are  two  ways  to  simulate  an  open  boundary.  One  may  either 
surround  the  domain  with  an  artificial  absorbing  medium,  so  that  outgoing  waves 
are  diffused  to  zero  before  they  return  to  the  computational  domain,  or  one  may  use 
a  differential  operator  to  prescribe  the  wave  behavior  at  the  boundary,  so  that  only 
outgoing  waves  are  permitted. 

Research  into  modeling  open  boundaries  has  been  active  since  the  late  1960s. 
Zienkiewicz  and  Newton  [121]  first  derived  the  Sommerfeld  radiation  condition  for 
outgoing  wave  propagation  in  1969.  In  hindsight,  this  differential  operator  is  surpris¬ 
ingly  obvious  and  easy  to  derive.  Beginning  with  the  known  general  solution  to  the 
one-dimensional  scalar  wave  equation 

dttu  =  c^dxxU  ,  (III.l) 
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we  have 


u{x,  t)  =  F[x  —  ct)  +  G{x  +  ct)  (III-2) 

for  some  functions  F  and  G  [105].  If  F  denotes  the  outgoing  waves,  and  G  the 
incoming,  then  we  insist  on  G  =  0  on  the  boundary.  Differentiating  u  with  respect 
to  X  and  t,  this  gives  us 

dtu  =  -cF'  ,  d^u  =  F'  .  (III.3) 

From  here  we  easily  get  the  Sommerfeld  boundary  condition: 

dtU  +  cdxU  =  0  . 

This  boundary  condition  can  be  interpreted  in  two  ways.  The  characteristic-based 
interpretation  uses  this  condition  as  prescribing  the  Riemann  invariant  of  the  solution 
(see  Ch.  8  of  [24]).  The  wave-based  interpretation  describes  it  as  requiring  waves  on 
the  boundary  to  satisfy  the  one-way  advection  equation.  Several  early  NRBCs  were 
developed  in  the  1970s  using  these  two  interpretations.  Wurtele  et  ah  [119]  used 
the  characteristic  method,  while  Pearson  [93]  and  Orlanski  [92]  took  the  wave-based 
approach.  Engquist  and  Majda  [25,  26]  extended  the  wave-based  method,  dehning 
a  pseudo-differential  operator  solution  to  the  2-D  wave  equation  and  deriving  Fade 
approximations  thereto  in  a  sequence  of  ever-more-accurate  boundary  conditions. 
Smith  [101]  took  a  simplistic,  albeit  computationally  intensive,  approach:  Apply  a 
Dirichlet  boundary  condition  to  one  solution,  apply  a  Neumann  boundary  condition 
to  another,  and  then  add  the  two  solutions,  cancelling  the  two  reflections  and  leaving 
only  the  non-reflecting  solution  within  the  accuracy  of  the  Neumann  operator.  One 
of  the  hrst  absorbing  layer  methods  was  also  published  around  this  time  by  Davies 
[17]  for  the  linearized  Euler  equations  in  a  nested  environment,  using  a  “relaxation” 
function  near  the  boundary  to  match  the  small-scale  interior  scheme  with  the  large- 
scale  outer  model.  Later  NRBCs  built  on  these  methods  or  developed  new  approaches. 

In  the  1980s,  several  new  NRBC  techniques  appeared.  Bayliss  and  Turkel 
[9,  10]  developed  a  sequence  of  increasing-order  boundary  conditions  based  on  an- 
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nihilating  the  first  terms  of  the  asymptotic  expansion  of  the  scalar  wave  eqnation. 
Miller  and  Thorpe  [87]  proposed  a  modihed  Orlanski  scheme  based  on  alternative 
time-stepping  methods.  Perm  and  Gnstafsson  [27]  nsed  Fonrier  transformations  to 
devise  a  downstream  bonndary  condition  for  the  steady-state  linearized  Enler  eqna- 
tions.  Klemp  and  Dnrran  [79]  developed  a  “sponge  layer”  to  absorb  ontgoing  gravity 
waves,  applied  to  the  linear  Bonssinesq  eqnations,  later  applied  to  the  non-linear  Enler 
eqnations  by  Giraldo  and  Restelli  [32],  Davies  [18]  compared  varions  techniqnes,  in- 
clnding  “diffnsion  zones,”  relaxation  fnnctions,  and  radiation  bonndary  conditions. 
Raymond  and  Kno  [95]  modihed  the  Sommerfeld  condition  to  consider  tangential  as 
well  as  normal  derivatives  in  mnlti-dimensional  hows.  Ting  and  Miksis  [110]  pro¬ 
posed  nsing  Kirchhoh’s  formnla  to  determine  the  bonndary  valnes  of  waves  exterior 
to  a  scatterer.  Trefethen  and  Halpern  [111]  analyzed  the  Engqnist-Majda  method; 
considering  diherent  approaches  to  approximating  the  psendo-diherential  operator, 
they  demonstrated  its  well-posedness  and  proved  Engqnist’s  and  Majda’s  proposition 
concerning  which  classes  of  Fade  approximations  are  well-posed.  Higdon  [62,  64]  de¬ 
veloped  a  seqnence  of  increasing-order  bonndary  conditions  based  on  a  prodnct  of 
Sommerfeld  conditions  at  varions  incidence  angles.  Keller  and  Givoli  [78]  developed 
the  “Dirichlet-to-Nenmann”  (DtN)  mapping,  an  NRBG  method  which  converts  the 
Dirichlet  condition  at  inhnity  to  a  Nenmann  condition  at  the  compntational  domain 
bonndary;  they  then  nsed  this  DtN  mapping  in  a  hnite  element  solntion  of  Laplace’s 
eqnation  on  an  inhnite  domain  [38].  See  also  Givoli’s  1991  paper  [33]  reviewing  the 
then-cnrrent  state  of  the  art. 

The  1990s  saw  an  explosion  in  NRBG  development.  Kroner  [82]  adapted  the 
Engqnist-Majda  scheme  to  the  2-D  linearized  Enler  eqnations,  nsing  Fonrier  transfor¬ 
mations  rather  than  psendo-diherential  operators  to  dehne  the  bonndary  condition; 
Giles  [30]  performed  a  similar  techniqne,  also  nsing  Fonrier  transformations  to  de¬ 
rive  an  NRBG  for  the  2-D  linearized  Enler  eqnations.  G.  Kreiss  [80]  nsed  a  simple 
Dirichlet  condition  at  the  downstream  end  for  the  pressnre  term  of  the  linearized 
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Euler  equations,  and  he  showed  that  the  resulting  error  decreases  exponentially  up¬ 
stream  of  the  open  boundary.  Collino  [15]  extended  the  Engquist-Majda  scheme  to 
open  domains  (requiring  consideration  of  corner  conditions)  and  to  other  wave-like 
equations.  Tam  and  Webb  [106]  used  the  asymptotic  expansion  to  dehne  radiation 
boundary  conditions  compatible  with  their  dispersion-relation-preserving  hnite  dif¬ 
ference  scheme  for  the  linearized  Euler  equations.  Higdon  continued  his  sequence 
of  NRBC  papers  [65,  66,  67,  68],  culminating  in  a  robust  NRBC  sequence  for  the 
dispersive  (Klein-Gordon)  wave  equation.  Grote  and  Keller  [47,  48,  49]  extended  the 
DtN  technique  to  spherical  waves  and  the  Helmholtz  equation,  in  hnite  difference  and 
hnite  element  methods;  Thompson  and  Huan  [109]  later  modihed  this  formulation  to 
implement  a  hnite  element  solution  of  the  spherical  wave  equation.  Ren  [96]  used  a 
2-D  Sommerfeld-like  boundary  condition, 

{dt  + c{axd^  + aydy))r]  =  0  (HI. 5) 

=  1  , 

to  dehne  an  open  boundary  for  the  2-D  Boussinesq  equations.  Hagstrom  and  Hariha- 
ran  [54]  presented  an  NRBG  for  the  2-D  and  3-D  wave  equation  in  polar/spherical  co¬ 
ordinates,  using  auxiliary  variables  to  remove  the  high-order  normal  derivative  terms; 
Huan  and  Thompson  [75]  later  modihed  this  scheme  by  using  a  spherical  harmonic- 
based  formulation.  Safjan  [98]  also  took  the  auxiliary  variable  approach,  applying 
them  to  high-order  Fade  approximations  to  the  pseudo-differential  operator  of  the 
scalar  wave  equation.  Jensen  [77]  compared  several  techniques,  including  NRBGs 
and  sponge  layers,  for  modeling  open  boundaries  in  a  stratihed  ocean  model. 

The  1990s  also  saw  the  development  of  the  Perfectly  Matched  Layer  (PML), 
hrst  dehned  by  Berenger  [11]  for  the  2-D  Maxwell  equations.  This  absorbing  layer 
method  surrounds  the  computational  domain  with  a  dispersive  medium,  dehned  in 
such  a  way  that  incoming  waves  at  any  incidence  angle  pass  from  the  interior  to 
the  dispersive  layer  without  any  (theoretical)  rehection.  A  2007  paper  by  Skelton  et 
al.  [100]  claimed  to  hnd  over  1,000  references  to  Berenger’s  paper  in  the  literature. 
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This  technique  has  been  applied  to  the  linearized  Euler  equations  by  Abarbanel  et 
ah  [1],  Goodrich  and  Hagstrom  [45],  Hu  [71,  72,  73],  and  Natal  [88];  to  the  linearized 
shallow-water  equations  by  Navon  et  ah  [89];  to  Maxwell’s  equations  using  a  second- 
order  discretization  scheme  by  Sjogreen  and  Petersson  [99];  to  the  linearized  and 
non-linear  wave  equation  by  Appelo  and  G.  Kreiss  [7]  (following  Appelo  et  ah’s  well- 
posedness  and  stability  theory  in  [6]);  to  elastic  waveguides  by  Skelton  et  ah  [100]; 
to  the  time-harmonic  wave  equation  by  Bermudez  et  ah  [12];  and  to  the  non-linear 
Euler  and  Navier-Stokes  equations  by  Hu  et  ah  [74], 

NRBG  development  continued  in  this  decade  as  well,  with  new  techniques  and 
new  applications  of  old  techniques.  Oliveira  [91]  combined  a  Sommerfeld  condition 
with  an  absorbing  layer  to  develop  an  NRBG  for  the  transient  mild-slope  equation: 

V  [CCgVr])  +  {eCCg  -  a;2)  77  -  dui]  =  0  (HI.6) 

(see  Eqn.  (9)  of  [91]).  Grote  and  Keller  [50]  developed  an  exact  NRBG  for  the  3-D 
elastic  wave  equation  with  a  spherical  open  boundary,  based  on  annihilating  the  hrst 
terms  of  the  wave’s  spherical  harmonics.  Alpert  et  ah  [4,  5]  developed  two  NRBGs 
for  the  scalar  wave  equation,  one  using  Hankel  functions  and  Laplace  transforms,  and 
one  using  Fourier-Laplace  transformations.  Lie  [83]  used  Fourier-Laplace  transfor¬ 
mations  to  derive  an  NRBG  for  the  shallow-water  equations.  Golonius  and  Ran  [16] 
developed  an  absorbing  buffer  technique  for  conservation  law-based  systems  by  using 
Fourier  transformations  to  filter  the  outgoing  flow  disturbances.  McDonald  [85,  86] 
derived  a  characteristic-based  NRBG  for  the  shallow- water  equations  and  compared 
the  performance  of  NRBGs  and  “relaxation  zone”  boundaries  in  nested-model  envi¬ 
ronments.  Hagstrom  and  Nordstrom  [56]  analyzed  the  use  of  extrapolation  boundary 
values  in  solving  the  steady-state  linearized  Euler  equations,  showing  the  relation¬ 
ship  between  the  position  of  the  artihcial  far-held  boundary  and  the  error  norm  of 
the  discrete  solution.  Blayo  and  Debreu  [13]  considered  a  characteristic  variable  ap¬ 
proach  to  NRBGs  in  hrst-order  systems  for  ocean  and  atmospheric  models.  Ghang 
et  ah  [14]  used  a  Space-Time  Gonservation  Element  and  Solution  Element  (GE/SE) 
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method  to  solve  the  1-D  Euler  equations,  which  effectively  handled  shockwaves  inside 
the  domain  (although  the  results  were  good  within  the  domain,  the  numerical  results 
published  showed  poorer  performance  closer  to  the  boundary).  Guddati  and  him  [51] 
used  continued  fractions  (rather  than  the  Engquist-Majda  Fade  approximations)  to 
approximate  the  pseudo-differential  operator,  and  they  devised  a  formulation  that 
could  be  applied  to  any  convex  polygon  boundary  (rather  than  the  usual  straight-line 
boundaries).  Zahid  and  Guddati  [120]  incorporated  PML-like  “padding  elements” 
into  a  continued  fraction  NRBG  for  dispersive  waves.  Atassi  and  Galan  [8]  used 
Fourier-Bessel  modes  to  derive  an  NRBG  for  the  non-linear  Euler  equations  in  an 
annular  duct.  Song  and  Bazyar  [102]  derived  an  NRBG  for  hnite  element  frequency- 
domain  wave  analysis  based  on  Fade  approximations. 

The  2000s  also  saw  a  revived  interest  in  the  Higdon  NRBG  sequence.  Givoli 
and  Neta  created  an  algorithm  to  compute  high-order  hnite  difference  derivatives 
automatically,  removing  the  algebraic  complexity  which  limited  the  original  Higdon 
sequence  to  third-order.  This  automation,  along  with  an  auxiliary  variable  method 
which  removed  the  high-order  normal  derivatives,  was  applied  to  the  Klein-Gordon 
equation  and  the  shallow-water  equations  in  a  sequence  of  papers  by  them  and  their 
students  and  colleagues  [39,  40,  41,  42,  43,  90,  113,  114,  115,  116].  Givoli  again 
published  a  review  of  current  NRBG  techniques  in  [35].  In  addition,  Hagstrom  and 
Warburton  [57]  developed  a  modihed  form  of  the  Givoli-Neta  auxiliary  variable  NRBG 
for  the  scalar  wave  equation.  This  new  method  was  expanded  and  analyzed  in  [37, 
53,  55,  84]. 

In  this  dissertation,  we  use  the  Givoli-Neta  automation  to  apply  the  Hig¬ 
don  NRBGs  to  the  linearized  Euler  equations;  we  also  extend  the  Givoli-Neta  and 
Hagstrom- Warburton  auxiliary  variable  NRBGs  thereto.  Some  of  the  results  pre¬ 
sented  in  the  following  chapters  have  been  published  or  submitted  for  publication 
[19,  20,  21,  22].  The  remainder  of  this  chapter  is  devoted  to  the  scalar  form  of  the 
Higdon,  Givoli-Neta,  and  Hagstrom- Warburton  NRBG  sequences,  laying  the  ground- 
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work  for  their  application  to  the  linearized  Enler  eqnations. 


B.  NRBCS  FOR  SCALAR  EQUATIONS 
1.  Higdon 

One  of  the  simplest  non-reflecting  bonndary  conditions  (NRBCs)  is  the  Som- 
merfeld  radiation  condition: 

{dt  +  c,d,)u  =  Q,  (III.7) 

where  u  is  the  nnknown  solntion  to  onr  problem,  and  Cq  is  the  wave  propagation  speed 
in  the  positive  x  direction  (for  a  right-side  bonndary;  on  other  sides,  replace  with 
the  appropriate  normal  derivative).  In  essence,  this  bonndary  condition  says  that  the 
ontgoing  wave  at  the  bonndary  satishes  the  one-dimensional  advection  eqnation  with 
advection  speed  Cq.  This  bonndary  condition  is  most  easily  applied  to  the  standard 
wave  eqnation, 

duu  =  clV\  (III.8) 

It  can  also  apply  to  the  linearized  shallow  water  eqnations  (see  [41]), 

dtu  -  fv  =  -gd^r] 

dtv  +  fu  =  -gdyT]  (III. 9) 

dt-q  +  ho  {d^u  +  dyv)  =  0 

whose  snrface  height  component  rj  can  be  converted  (see  [113])  to  the  dispersive 
(Klein-Gordon)  wave  eqnation, 

dttV  =  clV'^rj  -  frj  ,  co  =  \fgho  (III.  10) 

The  difficnlty  in  (III. 9)  and  (III.  10)  is  that  there  is  more  than  one  wave  speed.  The 
same  problem  afflicts  the  standard  wave  eqnation  in  more  than  one  dimension,  as 
waves  travelling  in  directions  other  than  normal  to  the  bonndary  have  wave  speeds 
whose  normal  components  are  not  eqnal  to  cq.  Higdon  proposed  [68]  dehning  a 
bonndary  condition  as  a  prodnct  of  J  Sommerfeld-like  terms: 

j 

n(^i+cA)«=o  (III.  11) 
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Just  as  “a:?/  =  0”  is  a  consolidation  of  “a;  =  0  or  ?/  =  0” ,  this  product  consolidates,  into 
a  single  equation,  numerous  Sommerfeld  conditions  with  different  advection  speeds. 
Since  the  differential  operator  is  linear,  it  is  easy  to  show  that  if  any  one  of  the 
individual  Sommerfeld  conditions  is  satished,  then  the  consolidated  Higdon  condition 
is  also  satished. 

a.  Reflection  Coejficient 

Let  us  analyze  the  claim  that  this  boundary  condition  is  non-rehecting. 
Consider  a  wave-like  equation  with  a  solution  of  the  form 

u(a;,t)  (III.  12) 


This  solution  dehnes  a  wave  with  speed  Cx-  Applying  the  Higdon  boundary  condition 
to  this  equation,  we  claim  that 

n  =  0  (HI.  13) 

j=i 

If  Cx  does  not  exactly  equal  one  of  the  Cj,  this  statement  may  not  be  true.  To  make 
it  true,  we  must  add  a  rehected  wave,  thus: 

j 

n  +  Cjdx)  -f  =  0  (HI.  14) 

i=i 

The  magnitude  of  Rj  is  the  reflection  coefficient  associated  with  the  boundary  con¬ 
dition.  (If  (HI.  13)  is  true,  then  Rj  =  0.)  Higdon  [68]  claims  that  this  rehection 
coefficient  is  always  less  than  one,  and  that  it  decreases  as  J  increases.  We  present 
the  proof  here  in  more  detail. 


Lemma  III.l  The  magnitude  of  the  reflection  coefficient  Rj  of  the  Higdon  NRBC 
defined  in  (III.  11),  applied  to  a  wave  with  speed  Cx,  is 


J 


=  n 


C 


Proof.  We  prove  by  induction.  For  the  J  =  1  case,  we  have 


(HI.15) 


{dt  J-  cffix) 


0, 


(HI.16) 
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and  it  is  easy  to  show  that 


|i?i|  =  ^  .  (III.17) 

+  Cl 

Having  proven  there  exists  a  J  snch  that  this  lemma  is  trne,  we  consider  the  J  +  1 
case: 

J+i 

n  (dt  +  c,d^)  =  0  (111.18) 

i=i 

Since  the  differential  operators  are  linear,  we  can  write 

n  (dt  +  C,d^)  [(a*  +  =  0  .  (III.19) 

i=i 

Expanding  the  derivatives  inside  the  brackets  and  combining  terms,  we  get 

n  (dt  +  cjd,)  =  0  .  (III.20) 


(III.21) 


(III.22) 


Since  these  Cj  and  Cx  are  positive,  each  term  in  the  prodnct  is  less  than  one.  Thns,  as  J 
increases,  the  reffection  coefficient  gets  smaller.  Of  conrse,  well-chosen  Cj’s  will  rednce 
the  reffection  coefficient  more  rapidly,  bnt  even  poor  choices  (for  example,  failing  to 
consider  possible  incidence  angles  or  phase  speeds)  will  still  bring  improvements. 
Hence,  we  can  get  good  absorption  either  with  a  high  J  or  with  well-chosen  c/s. 
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b.  A  Simplifying  Assumption 

The  difficulty  in  using  (III.  11)  comes  from  the  rapidly  increasing  alge¬ 
braic  complexity  for  large  J.  At  J  =  2,  we  have 


^2 


,  ^  dtdx 


+  C2 


dtdx 


+  C1C2 


dx‘^ 


u  =  0 


(111.23) 


At  J  =  3,  it  becomes 


( 4-  r  r  r 

^  “I"  "I"  ^^dt^dx 


+ClC2g^  +  ClCsg^  C2C3g^  +  CiCsCs^) 


U 


y  =0 


(III.24) 


and  so  forth.  The  Givoli-Neta  algorithm  automates  the  hnite  difference  computation 
of  these  high-order  derivatives,  but  at  a  cost  of  0{3'^)  operations  per  time  step. 
However,  if  we  simply  make  all  cj  equal  to  a  single  wave  speed  c,  then  (III.  11)  can  be 
written  as 

{dt  +  cd^yu  =  0  (III.25) 


(see  [90]),  or,  expanding  the  polynomial, 

f  J  71  d'^  \ 

FxTTf— =  0  (111.26) 

We  will  show  in  the  next  section  that  this  simplihcation,  when  applied  to  a  hnite 
difference  formulation,  requires  only  O(J^)  operations  per  time  step. 

c.  Discretization  of  Higdon  NRBCs 

For  our  hnite  diherence  implementation  of  the  NRBCs,  we  use  hrst- 
order  diherences  for  each  derivative,  with  the  diherence  pointing  into  the  domain  and 
backward  in  time.  Higdon  [68]  demonstrated  that  this  discretization  scheme  is  stable 
when  used  in  conjunction  with  the  Klein-Gordon  equation.  For  an  open  boundary  on 
a  rectangular  domain,  we  have 
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0 


(111.27) 


Bottom:  o",  =  0, 

where  I  denotes  the  identity  operator;  S~  and  S~^  backward  and  forward  shifts  in  the 
snbscript  variable,  e.g., 


0+^77  _ 


s; 


a. 


=  cr. 


n— 1. 
’ 


St  onr  grid  spacing  in  time;  Sx  and  Sy  onr  grid  spacing  in  x  and  y,  respectively;  crfj 
the  valne  of  a  generic  state  variable  a  at  grid  point  {i,j)  at  time  step  n;  and  snbscripts 
N,  W,  E,  S  denoting  the  north,  west,  east,  or  sonth  bonndaries,  respectively.  If  we 
mnltiply  each  term  by  St,  clearing  it  ont  of  the  denominator,  and  then  gronp  terms 
by  shift  operator,  we  get 


Top: 

(^ttyl  +  bS^  +  CySy  ^ 

=  0 

Left: 

{ojxI  +  bSj-  +  CxS'^^  a^j 

=  0 

Right: 

(ttxl  +  bS^  +  )  a%  j 

=  0 

(III.28) 

Bottom: 

{ttyl  +  bS-  +  CyS+y  a^s 

=  0, 

where 


ttr  =  1  —  Cr 


Qjy  1  Cy 


b  = 


('X  Cq 


Cy 


a 

Sx 
St 


Expanding  these  operators  as  polynomials,  we  have 

Top:  ( E  E  cist's;-’ 


cr. 


i,N 


0 
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Left: 

Right: 

Bottom: 


j  J-0 


EE 

3=0  7=0 
(  J  J-0 

EE 


J! 


J! 


^^0  ^0 
J  J-0 

EE 


E,j 


J! 


fr«c2S,-'’sr  < 


0 

0 


(111.29) 


where  a  =  J  —  (3  —  'y. 

Note.  From  this  point  forward  in  the  dissertation,  we  will  only  show 
the  NRBC  formnla  for  one  side  rather  than  all  fonr.  On  the  other  three  sides,  the 
appropriate  changes  shonld  be  made  for  the  correct  normal  derivative. 


2.  Givoli-Neta 

The  greatest  problem  afflicting  these  high-order  NRBCs  is  the  presence  of 
high-order  spatial  and  temporal  derivatives.  With  the  spatial  derivatives,  the  NRBC 
algorithm  mnst  look  deep  into  the  domain,  and  an  incoming  wave  can  begin  to  affect 
the  bonndary  long  before  it  actnally  reaches  the  bonndary.  High  temporal  derivatives 
reqnire  a  long  “history”  of  past  valnes,  increasing  memory  reqnirements.  In  both 
cases,  and  with  the  high-order  mixed  derivatives,  increasing  the  nnmber  of  terms  in 
the  NRBC  calcnlation  increases  the  danger  of  ronnd-off  errors  corrnpting  the  solntion 
and  destabilizing  the  system. 

In  addition,  the  NRBC  order  is  inherently  limited  by  the  size  of  the  domain. 
This  is  trne  in  a  hnite  difference  setting,  bnt  it  is  even  more  restrictive  in  dealing  with 
hnite  elements.  In  a  hnite  element  scheme,  the  spatial  derivatives  generally  limit  the 
NRBC  order  to  that  of  the  element  polynomials,  unless  one  is  willing  to  nse  derivative 
approximations  which  are  not  local  to  a  single  element. 

The  way  ont  of  this  dilemma  is  to  introdnce  auxiliary  variables  [34]  which 
reqnire  only  low-order  derivative  calcnlations.  In  this  section  and  the  one  following, 
we  consider  two  methods  for  implementing  these  anxiliary  variable  techniqnes. 

The  Givoli-Neta  NRBC  was  hrst  proposed  in  [39,  42,  43]  for  the  Klein-Gordon 
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equation  in  finite  difference  and  finite  element  schemes.  (In  this  dissertation,  the 
term  “Givoli-Neta  NRBC”  refers  to  the  auxiliary  variable  implementation,  not  their 
automated  high-order  Higdon  NRBC.)  This  NRBC  follows  directly  from  the  Higdon 
NRBC.  Where  the  Higdon  NRBC  uses  one  high-order  equation  on  the  boundary,  the 
Givoli-Neta  NRBC  of  order  J  uses  a  system  of  J  low-order  equations,  dehned  thus: 

<fj+i  =  ‘fj  (HI. 30) 

(po  =  u 

(pj  =  0  . 

Direct  substitution  shows  that  (HI. 30)  is  equivalent  to  (HI.  11).  The  above  dehnition 
applies  to  an  NRBC  on  the  right  side  of  the  domain;  on  other  sides,  replace  dx  with 
the  appropriate  normal  derivative. 

To  illustrate  the  utility  of  this  formulation,  let  us  apply  the  NRBC  to  the 
Klein-Gordon  equation  (HI.  10).  As  this  work  has  already  been  done  by  Givoli  and 
Neta  [39,  42,  43],  we  will  not  repeat  the  derivations  here. 

It  is  easy  to  show  that  the  auxiliary  variables  also  satisfy  the  Klein-Gordon 
equation,  that  is, 

dtt^j  =  CoVVi  -  /Vi  ,  (HI.31) 

for  all  J  G  1 ...  J  —  1.  Armed  with  that  fact,  we  can  then  combine  (HI. 30)  and  (HI.  10) 
to  remove  the  normal  derivative  terms  from  our  NRBC.  When  we  do  so,  our  resulting 
NRBC  is 

(  72  ~  3  )  ^tt'Pj-i  A  ( ~  I"  ~  )  ~  P^j+i  ~  ^yyP’j-i  ~  •  (HI. 32) 

\Co  C/  VC  C+i/  C 

This  equation  gives  us  a  tri-diagonal  system  of  equations  to  solve  for  our 
auxiliary  variables.  With  an  efficient  matrix  solver,  the  effort  required  is  0{J),  rather 
than  the  O(J^)  required  for  the  original  formulation,  and  it  is  0(J)  even  without 
setting  the  Cj  values  to  Cq.  Furthermore,  since  the  auxiliary  variables  are  dehned  solely 
using  temporal  and  tangential  derivatives,  we  can  restrict  them  to  the  boundary. 
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reducing  the  amount  of  memory  needed  to  store  their  values.  Finally,  since  the 
NRBC  uses  no  derivatives  beyond  second  order,  they  can  be  easily  applied  to  a  hnite 
element  system  to  any  order.  There  are  three  downsides,  however.  NRBCs  on  two 
adjacent  sides  require  consideration  of  corner  compatibility  conditions  (due  to  the 
presence  of  tangential  derivatives),  their  numeric  stability  is  less  robust  over  long  time- 
integrations  [37],  and  the  conversion  to  remove  the  normal  derivatives  introduces  an 
additional  source  of  discretization  error  which  produces  an  “error  floor”  that  cannot 
be  overcome  by  simply  increasing  J  [36]. 

3.  Hagstrom-Warburton 

The  Hagstrom-Warburton  NRBCs  were  hrst  presented  in  [57]  for  the  scalar 
wave  equation,  with  subsequent  extensions  and  analysis  in  [37,  53,  55].  This  NRBC 
scheme  also  uses  a  system  of  low-order  auxiliary  variable  equations  instead  of  a  single 
high-order  boundary  equation.  The  Hagstrom-Warburton  NRBC  of  order  J  is  given 
by 


dt^pi  =  aodtu  +  coda;U 
ajdtipj+i  -  cod^ipj+i  =  ajdtipj  +  cod^ifj 
^j+i  =  0  , 


(HI.33) 


where  the  parameters  aj  G  (0, 1]  are  chosen  by  the  user.  Hagstrom  and  Warburton 
show  that  this  dehnition,  applied  to  a  plane  wave  traveling  at  speed  cq  and  incidence 
angle  6,  results  in  a  reflection  coefficient  of 

Oo  —  cos 0  ( dj  —  cos 6\^ 

oo  +  cos  9  y  ttj  +  cos  9  J  ’ 


R  = 


(HI.34) 


which,  like  the  Higdon  scheme’s  reflection  coefficient,  is  a  product  of  J  terms,  each 
of  which  is  less  than  unity.  Hence,  as  J  increases,  the  reflection  coefficient  decreases. 
However,  with  the  squaring  of  each  term  in  the  product,  the  reffection  coefficient 
for  the  Hagstrom-Warburton  scheme  decreases  signihcantly  more  rapidly  than  the 
Higdon  reflection  coefficient. 
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No  physical  interpretation  of  the  NRBC  formulation  is  given.  One  way  to 
interpret  (III. 33b)  is  that  the  outgoing  characteristic  of  (pj  is  matched  by  the  incoming 
characteristic  of  pj+i.  However,  that  does  not  explain  the  form  of  (III. 33a).  It  could 
be  that  the  authors  wanted  an  NRBC  scheme  with  a  quadratically-decaying  reflection 
coefficient,  and  they  reverse-engineered  this  NRBC  formulation  to  get  one  that  works. 

As  with  the  Givoli-Neta  NRBCs,  the  auxiliary  variables  presented  here  also 
satisfy  the  Klein-Gordon  equation  (III. 31)  for  all  j  G  1 . . .  J.  We  can  use  that  fact  to 
remove  the  normal  derivative  term  from  (III. 33b),  replacing  that  equation  with 


-  l)  +  aj-i  (aj  -  l)  duPj+i 

—  (aj_i  -f-  ttj)  +  1)  duPj 


>  =  < 


(/Vi-l  -  cldyyPj-l) 

+  (%-l  +  %)  (/Vi  -  Cldyypj) 

+aj_i  (/Vi+1  -  cldyyPj+i) 

(III.35) 

(again,  see  [55]  for  details).  Once  again,  we  have  a  system  of  auxiliary  variables 
dehned  solely  on  the  boundary,  and  a  boundary  condition  consisting  entirely  of  low- 
order  derivatives.  As  with  the  Givoli-Neta  system,  this  system  also  requires  0{J) 
operations  to  solve. 


C.  NRBCS  FOR  FIRST-ORDER  SYSTEMS 

In  general,  the  three  NRBC  techniques  discussed  in  the  previous  section  have 
only  been  implemented  for  scalar  wave-like  equations,  not  for  hrst-order  systems. 
The  exception  is  the  Hagstrom-Warburton  NRBC,  where  one  section  of  [57]  briefly 
discusses  the  characteristic-based  implementation  for  a  symmetric  hrst-order  system 
of  the  form 

dtp  +  AdxP  +  Bdyp  =  0  .  (III. 36) 

However,  that  discussion  does  not  consider  corner  conditions  or  the  presence  of  un¬ 
differentiated  terms  on  the  right-hand  side.  In  Chapter  VI  we  develop  these  NRBCs 
in  more  depth. 
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While  the  Higdon  and  Givoli-Neta  NRBC  methods  have  been  developed  for 
the  shallow- water  eqnations  (e.g.,  [113]),  those  implementations  involve  converting 
the  system  to  a  scalar  Klein-Gordon  eqnation  for  the  snrface  height  state  variable; 
hence,  they  are  not  trnly  implemented  for  the  hrst-order  system.  In  the  next  two 
chapters,  we  will  develop  these  NRBG  methods  for  a  trne  hrst-order  system.  While 
the  Higdon  implementation  is  specihc  to  the  linearized  Enler  eqnations,  the  Givoli- 
Neta  implementation  (not  inclnding  gravitational  effects  (Sec.  V.G))  can  be  extended 
easily  to  any  hrst-order  system. 
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IV.  HIGDON  NRBCS  FOR  THE 
LINEARIZED  EULER  EQUATIONS 


A.  OVERVIEW 

This  chapter  develops  the  hnite  difference  implementation  of  the  high-order 
NRBCs  for  the  linearized  Enler  eqnations  in  two  dimensions.  We  begin  with  the 
simplest  possible  eqnation  set,  with  no  advection  or  forcing  terms.  After  demon¬ 
strating  the  validity  of  this  prototypical  implementation  in  Section  B,  we  proceed  to 
incorporate  the  effects  of  Coriolis  (Sec.  C),  gravity  (Sec.  D),  and  advection  (Sec.  E). 


B.  AN  INITIAL  PROTOTYPE 

We  begin  with  the  simplest  possible  implementation  of  the  linearized  2-D  Enler 
eqnations:  zero  advection,  no  Coriolis  or  gravitational  forces. 


9tP  +  Po  +  dyv)  =  0 

dtu  +  —dxP  =  0 
Po 

dtv  +  —dyp  =  0 

Po 

dtp  +  7Po  {dxU  +  dyv)  =  0 


(IV.l) 


The  variables  here  are  the  pertnrbation  variables  from  the  linearization  in  Sec.  II. C. 2; 
we  remove  the  asterisks  here  and  in  all  snbseqnent  eqnations  for  textnal  clarity.  In 
addition,  we  nse  po  and  po  instead  of  p  and  p  becanse  we  have  no  ^-dependency  in  these 
eqnations.  This  set  is  mostly  nseless  for  real-world  modeling,  especially  considering 
the  fact  that  p  has  no  impact  on  any  of  the  other  state  variables.  However,  it  will 
serve  as  an  initial  prototype  for  developing  this  implementation  (see  also  [19]). 

1.  Equivalence  of  (IV.l)  and  the  Scalar  Wave  Equation — 
Continuous  Case 

In  [68],  Higdon  proved  that  the  J-order  NRBC  is  compatible  and  stable  when 
applied  to  the  dispersive  wave  eqnation  (and  thns  the  standard  wave  eqnation  by 
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setting  /  =  0),  based  on  stability  criteria  developed  by  H.-O.  Kreiss  [81]  (with  a 
geometric  characteristic-based  interpretation  provided  by  Higdon  [63]).  Therefore,  if 
we  can  convert  these  simplihed  Enler  eqnations  to  the  standard  wave  eqnation,  we 
will  know  that  they  too  will  be  stable  with  this  NRBC  formnlation. 

Differentiate  (IV.lb,c,d)  with  respect  to  x,  y,  and  t,  respectively,  and  get 
(ignoring  the  eqnation  for  p) 

dxtu  dxxP 

Po 

dytV  =  -—dyyP  (lV.2) 

Po 

dttP  =  -IPO  {dxtu  +  dytv) 


Snbstitnting  the  hrst  and  second  of  these  new  eqnations  into  the  third  gives 

dttP  =  —  {dxxP  +  dyyp) , 

Po 

which  is  the  standard  wave  eqnation  for  p,  with  the  wave  speed  given  by 


(IV.3) 


Co  = 


npo 

po  ’ 


(IV.4) 


reflecting  the  valne  given  in  the  literatnre  [24,  76,  107,  108].  Therefore,  this  system  of 
eqnations,  combined  with  the  Higdon-like  NRBCs,  will  be  stable.  It  may  seem  a  bit 
shaky  to  make  this  claim,  since  we  have  only  proven  that  the  system  can  be  converted 
to  the  wave  eqnation  for  p,  not  for  the  other  variables.  Might  they  still  be  nnstable 
in  u  and  n?  The  answer  is  no.  Since  p  depends  on  u  and  n,  if  they  are  nnstable, 
then  they  will  destabilize  p.  Since  p  is  stable,  it  follows  that  u  and  v  mnst  also  be 
stable.  (The  stability  of  p  is  somewhat  irrelevant,  as  it  has  no  impact  on  the  other 
state  variables.) 

In  Sec.  VI.  1,  we  will  show  that  all  fonr  state  variables  satisfy  the  scalar  wave 
eqnation,  and  we  will  derive  its  exact  form.  For  now,  it  is  snfficient  to  note  that  we 
can  convert  this  eqnation  set  to  the  scalar  wave  eqnation  in  p. 
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2.  Equivalence  of  (IV.  1)  and  the  Scalar  Wave  Equation — 
Discrete  Case 

Using  stability  criteria  developed  by  Gustafsson  et.al.  [52],  Higdon  showed 
[68]  that  the  discretization  scheme  in  Sec.  III.B.l.c  is  compatible  with  the  standard 
second-order  discretization  scheme  for  the  standard  wave  eqnation  (solving  for 


u. 


n+l 


+  <f 


u. 


{Sty 


=  Cn 


2«m- 


+  <-i, 


(fe) 


+ 


u 


+Ki-y 


{Syy 


(IV.5) 


However,  this  scheme  is  for  solving  a  single  second-order  PDE,  not  a  system  of  hrst- 
order  PDEs,  as  we  have.  Even  thongh  we  have  shown  that  this  system  is  eqnivalent, 
in  one  of  the  state  variables,  to  the  standard  wave  eqnation,  creating  a  discretization 
scheme  which  is  also  consistent  and  stable  is  a  matter  of  some  delicacy. 

Givoli  and  Neta  described  a  similar  difficnlty  they  enconntered  in  developing 
a  stable  discretization  scheme  for  the  shallow  water  eqnations  [39].  After  nnmerons 
failed  efforts,  whose  resnlts  were  nnstable  for  any  J  >  2,  they  dehned  a  scheme  which 
was  eqnivalent  to  the  dispersive  wave  eqnation  discretization  scheme  proved  stable 
by  Higdon  in  [68].  We  nse  a  similar  approach  here.  Let  denote  a  second-order 
centered  difference  in  time,  with  similar  notation  for  difference  approximations  in  x 
and  y;  thns. 


A 


a 


a 


s*  -  s; 

25a 

e  {x,y,t}  , 


(IV.6) 


where  5a  is  the  step  size  in  a.  Use  the  following  discretization  for  (IV. 1): 


Atp 


Atti 


/\tv 


-po  {A^u  +  Ayv) 
A^p 


Po 

\P 

Po 

-7Po  {A^u  +  Ayv) 


(IV.7) 
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As  with  converting  (IV.  1)  to  (HI. 8),  apply  A^.  to  (IV. 7b),  Ay  to  (IV. 7c),  At  to  (IV. 7d), 
and  then  make  the  appropriate  snbstitntion.  This  gives  ns 


AtAtp  =  —  {A^A^p  +  AyAyp) 

Po 


(IV,8) 


If  we  expand  these  differences,  we  get 


Pit  -  ‘^Plj  +  Pli  _  7P0  ( Pl+2,j  -  ‘^Plj  +  Pl-2,j  Plj+2  -  “^Plj  +  Plj-2  \  ,  .X 

{2Sty  Po  I  {2Sxy  +  {2Syy  )  ^ 


which  is  the  scalar  wave  eqnation  on  a  double-sized  grid.  Hence,  to  make  onr  dis¬ 
cretization  scheme  fnlly  compatible,  we  nse  the  following  discretization  for  the  NRBC 
on  each  side  of  the  domain: 


i-sp\ 


(IV.  10) 


where  a  denotes  any  one  of  onr  state  variables  p,  u,  v,  p. 

A  reviewer  for  [21]  noted  that  this  NRBC  cannot  resolve  the  shortest  wave¬ 
lengths  resolvable  by  the  interior  scheme.  Snbseqnent  experiments  showed  no  new 
instabilities  generated  by  these  short  wavelengths. 

In  [19],  the  anthors  showed  that  a  certain  choice  of  one-sided  differencing  can 
also  yield  a  compatible  scheme,  one  which  is  compatible  with  the  scalar  wave  eqna¬ 
tion  on  the  same  size  grid.  However,  snbseqnent  experimentation  revealed  deeper 
flaws  in  that  approach.  First,  it  conld  not  be  extended  to  match  the  Klein-Gordon 
eqnation  when  Coriolis  forces  are  inclnded.  Second,  its  semi-implicit  natnre  became 
fnlly  implicit  nnder  advection.  (An  attempt  to  avoid  an  implicit  scheme,  based  on 
a  qnadrant-based  interior  discretization  method,  proved  to  be  nnstable  nnder  ad¬ 
vection.)  The  scheme  devised  here  avoids  these  flaws,  as  snbseqnent  sections  will 
show. 


3.  Numerical  Example — Semi-Infinite  Channel 

Consider  a  simple  nnmerical  example.  Dehne  a  sqnare  domain  10  km  on  each 
side,  with  walls  on  three  sides  and  open  on  the  top  (see  Fig.  6).  We  trnncate  the 
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domain  with  an  artificial  north  boundary  T  and  impose  the  non-reflecting  boundary 
condition  on  the  state  variables.  We  dehne  the  hard  wall  as 


dnP  =  0 

u-n  =  0  (IV.ll) 

dnP  =  0 

The  open  boundary  is  dehned  by  a  Higdon- type  NRBC  of  order  J.  Using  a  sea- level 


Figure  6.  A  semi-inhnite  channel  domain  truncated  by  an  artihcal  boundary  Ft 
[After  [40],  Fig.  lb,  p.  259] 

atmospheric  density  po  =  1-2^  and  pressure  po  =  1.01  x  10®^  [58],  our  initial 
condition  is  a  pressure  bubble  in  the  center  of  the  domain: 


Px,z  = 


To  1  + 


100 


d  <  r 


Po  :  otherwise 


pI,z 


=  < 


Po 


4^1"  :  d<r 


Pz 


Po  :  otherwise. 


(IV.12) 
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where 


and  {xc,  Vc)  denotes  the  center  of  the  domain.  The  initial  pertnrbation  of  p  is  dehned 
to  maintain  a  constant  potential  temperatnre  in  the  pressnre  pertnrbation  bnbble 
[24],  We  divide  the  domain  into  a  100  x  100  grid  and  rnn  the  simnlation  np  to  t  =  24 
s,  which  is  snfficient  time  for  the  wave  trongh  to  reach  the  corners.  We  also  rnn 
a  reference  solntion  on  a  10  x  20  km  domain  consisting  of  100  x  200  grid  points; 
this  reference  solntion  has  hard  walls  on  all  fonr  sides,  and  the  simnlation  dnration 
is  short  enongh  that  reflections  from  the  top  bonndary  will  not  re-enter  the  original 
domain.  By  nsing  a  compnted  reference  solntion  rather  than  an  analytic  solntion,  we 
can  attribnte  all  differences  to  the  NRBC  error  rather  than  the  discretization  scheme’s 
trnncation  error.  Onr  time  step  is  chosen  as  90%  of  the  CFL  limit  [64], 


2 


<1, 


(IV.13) 


where  Cq  is  the  wave  speed  \J^Po/ Po-  (We  have  observed  experimentally  [21]  that 
setting  5t  to  the  maximnm  allowed  by  the  CFL  limit  resnlts  in  rednced  effectiveness 
for  the  higher  J  NRBCs.)  We  nse  this  same  cq  for  the  NRBC  wave  speed.  Dehne  the 
error  norm,  for  each  state  variable  p,  to  be 


E^p  — 


(IV.  14) 


where  pj  is  the  state  variable  compnted  nsing  the  order  J  NRBC;  po  is  the  reference 


solntion  state  variable;  W  and  Ny  are  the  nnmber  of  points  in  the  x  and  y  directions, 
respectively;  and  we  normalize  the  error  norm  by  the  norm  of  the  reference  state 
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variable  (so  that  all  four  state  variables’  error  norms  are  approximately  the  same 
order  of  magnitude).  Using  the  perturbation  variables  dehned  in  (IV.  1),  the  interior 
discretization  scheme  (IV. 7),  and  the  Higdon  NRBC  discretization  (IV.  10a),  we  run 
the  simulation  using  different  values  of  J  for  the  NRBC  order.  Figs.  7-10  show 
the  perturbation  variables  p,  u,  u,  and  p,  respectively,  for  the  J  =  10  case.  Fig.  11 
contrasts  the  computed  solutions  and  error  deltas  for  v  between  the  J  =  1  and  J  =  10 
cases.  Table  I  shows  the  error  norms  for  each  state  variable  for  J  G  1 ...  10. 


Figure  7.  Plot  of  p  in  basic  system  (IV.  1)  with  J  =  10  in  a  semi-inhnite  channel. 
(TL)  Computed  solution.  (Center)  Reference  solution;  the  area  corresponding  to  the 
computed  solution  is  contained  below  the  horizontal  line.  (BL)  Reference  solution 
truncated  to  computed  solution  domain.  (BR)  Delta  between  reference  solution  and 
computed  solution,  with  error  norm  computed  by  (IV.  14). 
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xIO^ 


Figure  8.  Plot  of  u  in  basic  system  (IV.  1)  with  J  =  10  in  a  semi-infinite  channel. 
(TL)  Computed  solution.  (Center)  Reference  solution;  the  area  corresponding  to  the 
computed  solution  is  contained  below  the  horizontal  line.  (BL)  Reference  solution 
truncated  to  computed  solution  domain.  (BR)  Delta  between  reference  solution  and 
computed  solution,  with  error  norm  computed  by  (IV.  14). 


J 

Ep 

Eu 

Ev 

Ep 

1 

0.12361 

0.077449 

0.1674 

0.12361 

2 

0.039728 

0.020397 

0.047278 

0.039727 

3 

0.026079 

0.010879 

0.022051 

0.026078 

4 

0.022763 

0.0077222 

0.014192 

0.022763 

5 

0.021505 

0.0060029 

0.010178 

0.021505 

6 

0.020889 

0.0049066 

0.0078822 

0.020889 

7 

0.020551 

0.0041032 

0.0063328 

0.020551 

8 

0.020362 

0.0035019 

0.0052668 

0.020361 

9 

0.020249 

0.0030695 

0.0044665 

0.020249 

10 

0.020176 

0.0027737 

0.0038836 

0.020176 

Table  I.  Error  norms  for  basic  system  (IV. 1)  with  J  G  1 . . .  10  in  a  semi-inhnite 
channel 
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Figure  9.  Plot  of  v  in  basic  system  (IV.  1)  with  J  =  10  in  a  semi- infinite  channel. 
(TL)  Computed  solution.  (Center)  Reference  solution;  the  area  corresponding  to  the 
computed  solution  is  contained  below  the  horizontal  line.  (BL)  Reference  solution 
truncated  to  computed  solution  domain.  (BR)  Delta  between  reference  solution  and 
computed  solution,  with  error  norm  computed  by  (IV.  14). 
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Figure  10.  Plot  of  p  in  basic  system  (IV. 1)  with  J  =  10  in  a  semi-infinite  channel. 
(TL)  Computed  solution.  (Center)  Reference  solution;  the  area  corresponding  to  the 
computed  solution  is  contained  below  the  horizontal  line.  (BL)  Reference  solution 
truncated  to  computed  solution  domain.  (BR)  Delta  between  reference  solution  and 
computed  solution,  with  error  norm  computed  by  (IV.  14). 


0 

i 

Figure  11.  Comparison  of  v  in  basic  system  (IV.  1)  computed  with  J  =  1  and  J  =  10 
in  a  semi-inhnite  channel,  with  error  norms  computed  by  (IV.  14).  (TL)  Computed 
solution  for  J  =  1.  (TR)  Computed  solution  for  J  =  10.  (BL)  Delta  between  reference 
solution  and  J  =  1  computed  solution.  (BR)  Delta  between  reference  solution  and 
J  =  10  computed  solution. 
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4.  Numerical  Example — Open  Domain 

We  now  perform  the  same  simulation  in  an  open  domain  (see  Fig.  12),  using 
NRBCs  on  all  four  sides  for  the  state  variables.  Before  proceeding,  we  note  that  we 


Figure  12.  An  open  domain  12  truncated  by  artihcial  boundaries  F^,  F^,  Fi^  and  F b 
[After  [19],  Fig.  1,  p.  1] 

now  have  points  which  are  on  the  corners  of  two  open  boundaries,  and  we  pause  to 
address  this  potential  complication.  Looking  at  (IV.  10),  we  see  that  for  each  point 
on  an  open  boundary,  that  point  only  depends  on  itself  at  previous  times  and  on 
the  adjacent  interior  points.  The  boundary  points  are  independent  of  each  other. 
Furthermore,  due  to  the  discretization  (IV. 7)  of  the  interior  scheme,  there  are  no 
interior  points  which  depend  on  the  corner  points.  Fig.  13  illustrates  this  dependency 
for  the  top-right  corner  of  an  open  domain.  Notice  that  no  other  point  depends  on 
the  corner.  Hence,  our  boundary  condition  at  the  corner  is,  by  and  large,  irrelevant. 
For  our  implementation  here,  we  decree  that  the  corners  are  considered  part  of  the 
top/bottom  boundaries.  Our  reference  solution  is  a  30  x  30  km  domain  with  300x300 
grid  points,  situated  so  that  the  center  of  the  reference  domain  corresponds  to  the 
computational  domain.  When  we  run  this  simulation,  we  get  results  such  as  those  in 
Figs.  14,  15,  and  Table  II.  (From  here  on,  we  will  only  plot  the  results  for  one  state 
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Figure  13.  Interior  and  boundary  discretization  dependencies  for  points  near  the 
top-right  corner  of  an  open  domain.  The  black  arrows  show  the  interior  dependecies 
based  on  the  discretization  scheme  (IV.  7).  Blue  arrows  show  the  dependencies  of  the 
boundary  points  except  for  the  corner.  Green  arrows  show  the  dependency  if  the 
corner  is  considered  part  of  the  top;  red  arrows  show  the  dependency  if  the  corner  is 
considered  part  of  the  right. 

variable,  rather  than  all  four.)  Not  surprisingly,  given  the  symmetry  of  the  domain 
and  the  discretization  scheme,  the  errors  for  u  and  v  are  the  same. 


J 

E, 

Eu 

Ep 

1 

1.5544 

2.0918 

2.0918 

1.5558 

2 

0.4589 

0.61777 

0.61777 

0.45933 

3 

0.22861 

0.30055 

0.30055 

0.22882 

4 

0.15198 

0.19766 

0.19766 

0.15212 

5 

0.11326 

0.14588 

0.14588 

0.11336 

6 

0.089796 

0.11564 

0.11564 

0.08988 

7 

0.074067 

0.095798 

0.095798 

0.074136 

8 

0.06246 

0.082285 

0.082285 

0.062519 

9 

0.053427 

0.071617 

0.071617 

0.053476 

10 

0.046644 

0.062478 

0.062476 

0.046687 

Table  II.  Error  norms  for  basic  system  (IV.  1)  with  J  G  1 ...  10  in  an  open  domain 
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so  100  -so  0  so  100  ISO  200 


Error  norm-O  062470 


005 

0 
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Figure  14.  Plot  of  u  in  basic  system  (IV. 1)  with  J  =  10  in  an  open  domain.  (TL) 
Computed  solution.  (Right)  Reference  solution;  the  area  corresponding  to  the  com¬ 
puted  solution  is  contained  within  the  center  box.  (CL)  Reference  solution  truncated 
to  computed  solution  domain.  (BL)  Delta  between  reference  solution  and  computed 
solution,  with  error  norm  computed  by  (IV.  14). 


Figure  15.  Comparison  ofp  in  basic  system  (IV. 1)  computed  with  J  =  1  and  J  =  10 
in  an  open  domain,  with  error  norms  computed  by  (IV.  14).  (TL)  Computed  solution 
for  J  =  1.  (TR)  Computed  solution  for  J  =  10.  (BL)  Delta  between  reference 
solution  and  J  =  1  computed  solution.  (BR)  Delta  between  reference  solution  and 
J  =  10  computed  solution. 
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C.  CORIOLIS  FORCES  IN  THE  XV  PLANE 

Let  us  now  add  Coriolis  forces  to  our  equation  set  [21],  Leaving  the  advection 
terms  zero,  our  equation  set  is 


dtP  +  Po  (Xu  +  dyv)  =  0 
dtu  +  —dxP  =  fv 


Po 

dtv  H - dyp  =  -fu 

Po 

dtp  +  '^Po(d^u  +  dyv)  =  0  (IV.  15) 

This  set  is  somewhat  more  useful  for  atmospheric  modeling,  although  we  still  have  p 
decoupled  from  the  rest  of  the  system.  As  with  the  preceding  section,  we  now  wish 
to  show  that  (IV.  15)  is  equivalent  to  the  Klein-Gordon  equation  duP  =  CqV^P  —  /^p, 
which,  as  noted  previously,  Higdon  proved  is  stable  with  this  NRBC  formulation  [68]. 


1.  Equivalence  of  (IV.  15)  and  the  Klein-Gordon  Equation — 

Continuous  Case 

This  conversion  begins  the  same  as  in  the  previous  section.  Differentiate 
(IV.15d)  with  respect  to  t 


dttP  +  'fPoiXtU  +  dytv)  =  0 


(IV.16) 


Now  differentiate  (IV.  15b)  with  respect  to  x  and  (IV.  15c)  with  respect  to  y  and  add 


^XtU  T  OytV  (dxxP  T  ^yyP^  f  (^xU  OyU^ 

Po 


(IV.  17) 


Now  substitute  (IV.  17)  into  (IV.16) 


dttP  -  —  V^p  +  /7Po  (5„n  -  dyu)  =  0  . 
Po 


(IV.18) 


Differentiate  (IV.  15b)  with  respect  to  y  and  (IV.  15c)  with  respect  to  x  and  subtract 

( 


Oytu  dxtu  T 

Po 


9yxP  dxyP  f  (9yV  A  QxU^ 


(IV.19) 


V 


=0 
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Combine  terms  to  get 


/  {d^u  +  dyv)  =  -dt  {d^v  -  dyu)  .  (IV.20) 

Combine  (IV.15d)  and  (IV.20)  to  get 

fdtp  =  'ypodt  {da;V  -  dyu)  .  (IV.21) 

Integrate  (IV.21)  with  respect  to  time  to  get 

fip-  Po)  =  IPo  {dxV  -  dyu)  .  (IV.22) 

Finally,  substitute  (IV.22)  into  (IV.  18) 

dttP  -  —V^p  +  f{p-po)  =  0,  (IV.23) 

Po 

which  gives  us  the  Klein-Gordon  equation  for  the  pressure  perturbation  p  —  po,  again 
with  wave  speed  \j'ypo/Po-  As  an  aside,  we  notice  that  (IV.20)  involves  the  time 
derivative  of  the  curl  of  u  on  the  right-hand  side.  Thus,  if  /  =  0,  then  V  x  F  is 
constant. 


2.  Equivalence  of  (IV.  15)  and  the  Klein-Gordon  Equation- 
Discrete  Case 

We  would  like  to  hnd  a  discretization  scheme  which  is  equivalent  to  the  discrete 
Klein-Gordon  equation 


_  27/”  -I-  7;”t 

J 


.n—1 


U. 


St'^ 


=  Cn 


r+i.i  -  2<i + <-ij  ,  <i+i  “  + Kj-i 


SI 


+ 
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(IV.24) 

Given  the  similarity  to  the  discrete  scalar  wave  equation,  we  begin  by  using  the 
discretization  scheme  in  Sec.  B.2,  adding  the  Goriolis  terms  to  the  right-hand  side  to 
get 

^tP  =  -po  (AxM  +  ^yV) 

Atu  =  -^^  +  fv 


Po 

\P 

Po 


fu 


=  -7Po  (A^u  +  Ayv) 


Atv 

Atp 
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(IV.25) 


Apply  Aj.  to  (IV.25b),  to  (IV. 25c),  A^  to  (IV.25d),  and  then  make  the  appropriate 

snbstitntion.  This  gives  ns 


AtAtp  =  —  {A^A^p  +  AyAyp)  -  f'ypo  {A^v  -  Ayu) ,  (IV.26) 

Po 

Next,  apply  Ay  to  (IV. 25b)  and  A^,  to  (IV.25c),  then  snbtract  and  combine  terms  to 
get 

At  {AyU  -  A^v)  =  f  {A^u  +  Ayv)  (IV.27) 

We  snbstitnte  (IV.25d)  into  (IV.27)  to  get 

At  {AyU  -  A^v)  =  -^Atp  (IV.28) 

If  we  apply  At  to  (IV.26)  and  incorporate  (IV.28)  into  it,  we  get 


At  AtAtp - {A^A^p  +  AyAyp)  +  f  p  =0 

Po 


(IV.29) 


Thns,  the  qnantity  inside  the  brackets  is  constant  in  time.  Since  we  are  assnming 
a  closed  system  with  no  sonrce  fnnctions  (other  than  the  Coriolis  force  itself),  then 
this  qnantity  mnst  initially  be  zero  and  will  thns  always  be  zero.  Expanding  the 
differencing  terms  gives  ns 

^  in  (p^2,-^p?,+pi2, ,  , 

(25t)^  Po  \  {2SxY  (25|/)^  / 

(IV.30) 

As  with  eqnating  the  initial  basic  case  to  the  scalar  wave  eqnation,  we  now  have  a 
system  which  is  compatible  with  the  Klein-Gordon  eqnation  on  a  donble-size  grid.  So 
again,  onr  NRBC  will  be  given  by  (IV.  10). 


3.  Numerical  Examples 

For  this  example,  we  nse  the  same  domain  and  initial  conditions  as  with  the 
basic  system  in  Sec.  B.  Onr  Coriolis  force  /  is  based  on  an  angnlar  momemtnm 
O  =  7.292116  X  10“®  s“^  [112]  at  a  latitnde  (j)  of  30°  N,  thns. 


/  =  20sin(/)  =  7.292116  x  10“®  s“^ 
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J 

E, 

Eu 

Ey 

Ep 

1 

0.12361 

0.077448 

0.1674 

0.12361 

2 

0.039728 

0.020397 

0.047278 

0.039727 

3 

0.026079 

0.010879 

0.022051 

0.026078 

4 

0.022763 

0.0077222 

0.014192 

0.022763 

5 

0.021505 

0.0060029 

0.010178 

0.021505 

6 

0.020889 

0.0049066 

0.0078821 

0.020889 

7 

0.020551 

0.0041032 

0.0063328 

0.020551 

8 

0.020362 

0.0035019 

0.0052668 

0.020361 

9 

0.020249 

0.0030695 

0.0044665 

0.020249 

10 

0.020176 

0.0027737 

0.0038836 

0.020176 

Table  III.  Error  norms  for  Coriolis-influenced  system  (IV.  15)  with  J  G  1 ...  10  in  a 
semi-infinite  channel 


When  we  apply  this  force  to  our  domains,  we  get  results  which  are  virtually  indistin¬ 
guishable  from  those  in  the  preceding  section.  This  is  not  too  surprising,  since  the 
Coriolis  force  is  so  small  that  its  effects  are  unnoticeable  over  such  short  timeframes 
[28] .  Tables  III  and  IV  show  the  error  norms  for  all  four  state  variables  for  J  G  1 ...  10 
in  the  channel  and  open  domain,  respectively.  Suppose,  just  for  curiosity’s  sake,  we 
make  the  Earth  rotate  much  more  rapidly  and  increase  /  by  a  factor  of  a  thousand. 


J 

E, 

Eu 

Eu 

Ep 

1 

1.5544 

2.0917 

2.0917 

1.5558 

2 

0.4589 

0.61777 

0.61777 

0.45933 

3 

0.22861 

0.30055 

0.30054 

0.22882 

4 

0.15198 

0.19766 

0.19766 

0.15212 

5 

0.11326 

0.14588 

0.14588 

0.11336 

6 

0.089797 

0.11564 

0.11564 

0.08988 

7 

0.074067 

0.095798 

0.095797 

0.074136 

8 

0.062461 

0.082285 

0.082284 

0.062519 

9 

0.053427 

0.071617 

0.071617 

0.053476 

10 

0.046645 

0.062478 

0.062476 

0.046689 

Table  IV.  Error  norms  for  Coriolis-influenced  system  (IV.  15)  with  J  G  1 ...  10  in  an 
open  domain 
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J 

E, 

Eu 

Ey 

Ep 

1 

1.4093 

0.70281 

0.70253 

1.4033 

2 

0.4086 

0.20314 

0.20309 

0.40687 

3 

0.20321 

0.098226 

0.098178 

0.20235 

4 

0.13492 

0.064387 

0.064365 

0.13435 

5 

0.10056 

0.047464 

0.047443 

0.10013 

6 

0.079742 

0.037609 

0.037594 

0.079405 

7 

0.065817 

0.031154 

0.031144 

0.065539 

8 

0.055543 

0.026777 

0.026767 

0.055308 

9 

0.047521 

0.023331 

0.023321 

0.04732 

10 

0.041445 

0.020369 

0.020361 

0.04127 

Table  V.  Error  norms  with  J  G  1 ...  10  in  an  open  domain  nnder  artifically-large 
Coriolis  (IV.  15) 


to  /  =  0.07292116  s“^,  jnst  to  make  its  effects  noticeable.  If  we  do  that,  we  get  the 
error  norms  given  in  Table  V  for  the  open  domain.  A  comparison  of  the  resnlts  for 
small  and  large  J  is  exemplihed  in  Figs.  16  and  17.  If  we  look  at  the  u  and  v  plots 
side-by-side  (Fig.  18),  we  can  discern  the  clockwise  rotation  generated  by  the  Coriolis 
force. 
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X-velocity  perturbation,  J^l 


X-velocity  perturbation.  J=10 


X-vek>city  perturbation  error,  Jsi 


I 

0 

r 


X-vetocity  perturbation  error.  J«10 


I 

0 

i 


Figure  16.  Comparison  of  u  in  Coriolis-influenced  system  (IV.  15),  using  artificially- 
large  Coriolis,  computed  with  J  =  1  and  J  =  10  in  an  open  domain,  with  error  norms 
computed  by  (IV.  14).  (TL)  Computed  solution  for  J  =  1.  (TR)  Computed  solution 
for  J  =  10.  (BL)  Delta  between  reference  solution  and  J  =  1  computed  solution. 
(BR)  Delta  between  reference  solution  and  J  =  10  computed  solution. 


Y-velocity  perturbation.  J=1  Y-velocity  perturbation.  Jsio 


20  40  60  80  100  20  40  60  80  100 


Y-ve*ocity  perturbation  error,  J=1 


I 

I 


Y-vetocity  perturbation  error.  J=10 


0,05 

0 

-005 

•0.1 
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Figure  17.  Comparison  of  v  in  Coriolis-influenced  system  (IV.  15),  using  artihcially- 
large  Coriolis,  computed  with  J  =  1  and  J  =  10  in  an  open  domain,  with  error  norms 
computed  by  (IV.  14).  (TL)  Computed  solution  for  J  =  1.  (TR)  Computed  solution 
for  J  =  10.  (BL)  Delta  between  reference  solution  and  J  =  1  computed  solution. 
(BR)  Delta  between  reference  solution  and  J  =  10  computed  solution. 
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Figure  18.  Side-by-side  plot  of  u  and  u,  illustrating  the  rotation  generated  by  the 
Coriolis  acceleration.  The  superimposed  clockwise  arrows  highlight  the  combined 
result  of  the  u  and  v  values. 
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D.  GRAVITATIONAL  FORCES  IN  THE  AZ  PLANE 


Now  we  shift  from  the  xy  plane  to  the  xz  plane  and  consider  what  happens 
when  we  add  the  effects  of  gravity  to  onr  system.  Still  keeping  the  advection  terms 
eqnal  to  zero,  we  now  have  the  eqnation  set 


dtp  +  p  {d^u  +  d^w)  =  -p'w 

dtu  +  -d^p  =  0 

P 

dtw  +  -d^p  =  --p 

P  P 

dtP  +  7P  (Am  +  d^w)  =  gpw 


(IV.31) 


where  we  now  mnst  nse  p  and  p  for  onr  ^-dependent  reference  states.  We  nse  p'  to 
denote  the  derivative  of  p,  which  depends  only  on  Note  that  p  hnally  makes  an 
impact  on  the  other  state  variables. 

1.  Defining  the  Reference  State  for  Density  and  Pres¬ 
sure 

Eq.  (11.59)  sets  a  compatibility  condition  between  onr  reference  states  for  den¬ 
sity  and  pressnre.  This  restriction  is  a  fairly  loose  one,  however,  reqniring  ns  to  look 
to  other  sonrces  for  possible  fnnctions  which  satisfy  this  condition.  Since  these  eqna- 
tions  are  derived  from  the  ideal  gas  law,  we  look  to  the  literatnre  for  atmospheric 
models  on  which  we  can  base  these  initial  conditions. 

Althongh  several  atmospheric  models  exist  [24,  112],  for  ease  of  differentiation 
we  will  nse  an  exponentially-decaying  model 


—  — OlZ 

p  =  Poe 


(IV.32) 


where  a  is  a  scaling  height  needed  to  match  the  snrface  {z  =  0)  pressnre  and  density 
valnes.  Applying  (11.59)  to  this  model,  we  get 


P 


Pqo: 

- e 

9 


(IV.33) 


57 


2.  Deriving  a  Wave-Like  Equation  from  (IV. 31) 

As  with  the  basic  system  and  the  system  subject  to  Coriolis  forces,  we  now 
attempt  to  convert  the  gravity- affected  system  to  a  wave-like  equation.  First,  we 
differentiate  (IV.31d)  with  respect  to  time  to  get 

dttP  =  gpdtw  -  7p  {d^tu  +  d^tw)  (IV. 34) 


Differentiating  (IV. 31b)  with  respect  to  x  and  (IV. 31c)  with  respect  to  ^  (remember¬ 
ing  that  p  depends  on  z)  gives  us 


d^tu  =  --d^^p 
P 


dztw  =  -g 


'pd;,p-p'p\  pdzzp-p'd:,p 


P^ 


Substituting  these  results  into  (IV. 34)  gives  us 


P^ 


p' 


dttP  =  gpdtW  +  —  dz:xP  +  gdzp - —p  -  —d^p 

P  \  P  P  / 

,  -o  ,  (  fi  9P'  p'  a  \ 

=  —V  P  -f  gpdtW  +  —  gd^p - —p  -  —d^p 

P  P  \  P  P  ) 

Using  (IV. 31c)  to  remove  the  pdtw  term  from  (IV.37)  gives  us 

o  2  a  ^  ^P  (  a  9P'  P'  r.  \ 

OttP  =  —V  p-  g  p-  go^p  +  —  gd^p - —p  -  —d^p 


P 


P 


P 


P 


IP  2  (  ,  1PP'\  ,  IPf.  (  .  1PP'\  ^ 

=  yV  U  + ^  I  P  +  U  + ^  I 

If  we  solve  (IV. 31a)  and  (IV.31d)  for  dxU  +  dzW,  we  get 

gpw  -  dtp  _  -p'w  -  dtp 
IP  P 

which  we  can  solve  for  w  and  differentiate  with  respect  to  time  to  get 

pdttP  -  ipdttp 


dtw  = 


gp^  +  7pp' 


,  P  7^  0 


Equating  this  with  dtW  in  (IV.31c)  and  solving  for  p  gives 

P  ( pdttP  -  ipdttp  1  \ 

P  =  -  - -2  ,  — , - -dzP 

g  \  gp  +  ipp  p  J 


(IV.35) 

(IV.36) 

(IV.37) 

(IV.38) 

(IV.39) 

(IV.40) 

(IV.41) 
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Substituting  this  result  into  (IV. 38)  gives  us 


dttP  = 


f  W  -p{a  +  ^-f) 

[  +9fd,p-{g  +  af)d,p 

giPpdttp  -  gp'^dttP  7PP'  ipdttp  -  pdttP 


IP  ^2 
=  —  V  p  — 

p 


gp"^  +  7PP' 


p  7P^  +  7pp' 


+  p^a,p(IV.42) 


after  combining  and  cancelling  terms.  Now  let  us  use  an  exponential  atmospheric 
model.  Combining  (IV. 32)  and  (11.59)  gives  us 


P  = 


P  = 


pa 

g 


pa 


Substituting  these  values  into  (IV.42)  gives  us 


dttP  =  —V  p- 
P 

IP 


7P^2„  -  g$9uP  7y  ( idttp  -  ^dup'' 


Cp'  Qp' 

g^-iT 


g 


+ 


IP. 


\  3f-Y- 


V  p  -  '^-dttp  +  dttP  +  g—dzP 

pa  p 

Noting  that  p/p  =  g/a,  this  simplihes  to 


+  g—o^p 
P 


(IV.43) 


dttp  =  V^p  +  gd^p 


(IV.44) 


We  almost  have  a  wave-like  equation.  The  only  thing  missing  is  having  the  time 
derivative  and  the  Laplace  operator  in  the  same  variable.  Look  again  at  (IV.41).  If 
we  solve  it  for  duP  and  use  (IV.32)  to  cancel  and  combine  terms,  we  have 

dttp  =  —dttP  +  — — —  {gp  +  d^p)  (IV.45) 

IP  7 

Substituting  this  result  into  (IV.44)  gives  us 

dttP  =  ?  (V^p  -h  gd^p)  +  p  (1  -  7)  pdtw  ,  (IV.46) 

P  ^  ' 

where  we  use  (IV. 31c)  to  make  the  substitution 


pdtW  =  -gp  -  d^p  , 
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and  we  keep  the  usual  ^  as  the  wave  speed  term  rather  than  simplify  it.  If  we  try  to 
remove  the  p  and  w  terms,  we  go  around  in  circles,  never  quite  replacing  them  all  in 
terms  ofp.  Hence,  this  is  the  closest  we  can  come  to  deriving  a  wave-like  equation  that 
is  equivalent  to  (IV.31).  Our  wave  equation  for  p  is  altered  by  the  vertical  derivative 
of  p  and  the  time  derivative  of  w. 

Note,  however,  that  if  we  set  7  =  1,  differentiate  (IV.46)  with  respect  to  time, 
and  make  the  appropriate  substitutions,  we  can  get 

dt  (dttp  -  -  PdztP  =  0  ,  (IV.47) 

which  more  nearly  resembles  a  wave  equation,  only  for  the  time  derivative  of  the 
pressure,  rather  than  for  the  pressure  itself.  However,  due  to  the  dehnitions  of  the 
atmospheric  constants,  7  =  1  means  R  =  0,  which  causes  problems  with  the  ideal 
gas  law  p  =  pRT  on  which  we  base  our  model.  Hence,  this  contrivance  is  purely 
academic. 

3.  Numerical  Examples 

For  our  hrst  example,  we  consider  an  open  boundary  on  the  sides  of  the  do¬ 
main,  so  that  the  normal  derivative  is  perpendicular  to  the  force  of  gravity  (see 
Fig.  19).  We  discretize  (IV.31)  according  to  the  same  second-order  centered-difference 
scheme  in  Sec.  B,  which  gives  us 


\p  = 

-p{A^u 

<1 

+ 

—  p'w 

Atu  = 

A^p 

P 

Atw  = 

A^p 

9 

-P 

P 

P 

Atp  = 

-IP  ( 

<1 

+ 

s 

1  +  gpw 

where  we  use  z,  w,  p,  and  p  instead  of  y,  v,  po,  and  po,  respectively,  to  emphasize 
that  our  concern  is  now  with  a  vertical,  rather  than  horizontal,  plane.  As  before,  we 
use  the  same  domain  size  and  initial  conditions.  However,  we  now  have  hard  walls 
on  the  top  and  bottom  and  NRBCs  on  the  left  and  right,  and  the  reference  solution 
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£2 


Figure  19.  An  infinite  channel  domain  truncated  by  artifical  boundaries  and 
Tr  [After  [22],  Fig.  1,  p.  2] 

is  a  30  X  10  km  channel,  with  the  domain  of  interest  in  the  center.  The  error  norms 
are  given  in  Table  VI,  with  a  comparison  of  pressure  plots  given  by  Fig.  20. 

Having  successfully  implemented  gravity  perpendicular  to  the  open  bound¬ 
aries’  normal  vector,  we  now  implement  it  parallel  to  the  open  boundary’s  normal 
vector.  Using  again  the  same  initial  conditions,  our  domain  is  now  a  “bucket,”  i.e., 
a  semi-inhnite  channel  with  the  open  boundary  on  top  (Fig.  6).  Table  VII  shows  the 


J 

E, 

Eu 

Ejj 

Ep 

1 

0.23241 

0.52603 

0.12775 

0.23171 

2 

0.065394 

0.14839 

0.033298 

0.065251 

3 

0.03214 

0.069192 

0.017724 

0.032071 

4 

0.021193 

0.04438 

0.012525 

0.021144 

5 

0.015814 

0.031941 

0.0096649 

0.015776 

6 

0.012458 

0.024676 

0.007901 

0.012429 

7 

0.010184 

0.019914 

0.0066308 

0.01016 

8 

0.0086452 

0.016684 

0.0056677 

0.0086248 

9 

0.0076161 

0.01435 

0.0049472 

0.0075982 

10 

0.0068845 

0.012637 

0.0044394 

0.0068682 

Table  VI.  Error  norms  with  J  G  1 ...  10  in  a  horizontal  channel  with  gravitational 
forces  (IV. 31) 
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Pressure  perturbation.  J=1 


20  40  60  80  100 


Pressure  perturbation,  J=10 


20  40  60  80  100 


Pressure  perturbation  error.  J=1 


Pressure  perturbation  error.  J=10 


20  40  60  80  100 

Error  norm-0  0042545 


Figure  20.  Comparison  of  p  in  gravity-influenced  system  (IV.31)  computed  with  J  =  1 
and  J  =  10  in  an  infinite  horizontal  channel,  with  error  norms  computed  by  (IV.  14). 
(TL)  Computed  solution  for  J  =  1.  (TR)  Computed  solution  for  J  =  10.  (BL) 
Delta  between  reference  solution  and  J  =  1  computed  solution.  (BR)  Delta  between 
reference  solution  and  J  =  10  computed  solution. 


error  norms,  and  Fig.  21  compares  the  ^-velocity  perturbations  for  J  =  1  and  J  =  10. 

Next  we  combine  these  results  into  a  domain  open  on  three  sides,  with  a  hard 
wall  on  the  bottom  (see  Fig.  22).  Table  VIII  shows  the  error  norms,  and  Fig.  23 
compares  the  x-velocity  perturbations  for  J  =  1  and  J  =  10. 
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J 

Ep 

Eu 

Ey 

Ep 

1 

0.097404 

0.10278 

0.22495 

0.096936 

2 

0.034899 

0.027336 

0.063934 

0.034738 

3 

0.026286 

0.014095 

0.030089 

0.026195 

4 

0.024522 

0.0098734 

0.019451 

0.024443 

5 

0.023927 

0.00764 

0.013899 

0.023852 

6 

0.023651 

0.0062309 

0.010768 

0.023578 

7 

0.023506 

0.0052081 

0.0086306 

0.023433 

8 

0.023425 

0.0044326 

0.0071758 

0.023353 

9 

0.023378 

0.0038462 

0.0060809 

0.023307 

10 

0.023348 

0.00339 

0.0052913 

0.023276 

Table  VII.  Error  norms  with  J  G  1 ...  10  in  a  vertical  bncket  with  gravitational  forces 
(IV.31) 


J 

E, 

Eu 

Eu 

Ep 

1 

0.31226 

0.85386 

0.30744 

0.30865 

2 

0.088812 

0.24972 

0.089607 

0.087862 

3 

0.043421 

0.11951 

0.04392 

0.042972 

4 

0.028607 

0.077987 

0.029393 

0.028308 

5 

0.021264 

0.057307 

0.021724 

0.021042 

6 

0.016764 

0.04499 

0.017264 

0.016589 

7 

0.01378 

0.036859 

0.014216 

0.013637 

8 

0.011743 

0.031518 

0.012194 

0.01162 

9 

0.010292 

0.027596 

0.010665 

0.010185 

10 

0.0091922 

0.024287 

0.0094066 

0.0090963 

Table  VIII.  Error  norms  with  J  G  1 ...  10  in  an  open-air  domain  with  gravitational 
forces  (IV.31) 
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2-velocity  perturtaUon, 
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2-velocity  perturtation.  J=10 
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Figure  21.  Comparison  of  w  in  gravity-influenced  system  (IV. 31)  computed  with 
J  =  1  and  J  =  10  in  a  semi-infinite  vertical  channel,  with  error  norms  computed  by 
(IV. 14).  (TL)  Computed  solution  for  J  =  1.  (TR)  Computed  solution  for  J  =  10. 
(BL)  Delta  between  reference  solution  and  J  =  1  computed  solution.  (BR)  Delta 
between  reference  solution  and  J  =  10  computed  solution. 


n 


Figure  22.  A  half-plane  domain  D  truncated  by  artihcial  boundaries  Fr,  F^  and  Fr 
[After  [20],  Fig.  1,  p.  3] 
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Figure  23.  Comparison  of  u  in  gravity-influenced  system  (IV. 31)  computed  with 
J  =  1  and  J  =  10  in  an  open-air  domain,  with  error  norms  computed  by  (IV.  14). 
(TL)  Computed  solution  for  J  =  1.  (TR)  Computed  solution  for  J  =  10.  (BL) 
Delta  between  reference  solution  and  J  =  1  computed  solution.  (BR)  Delta  between 
reference  solution  and  J  =  10  computed  solution. 
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E.  ADVECTION 


Now  we  consider  the  case  of  non-zero  mean  flow.  Beginning  with  the  simple 
case,  we  have 


dtp  +  uod^p  +  vodyp  +  po  {d^u  +  dyv)  =  0 

dtu  4-  UodxU  +  VodyU  H - d^p  =  0 

Po 

dtv  +  u^dxV  +  vodyV  -\ - dyp  =  0 

Po 

dtp  +  uod^p  +  vodyP  -t-  7Po  +  dyv)  =  0 


(IV.49) 


Dehne  a  Lagrangian  derivative 


Ttt  —  dt  +  uodx  +  vody 
Dt 


If  we  nse  this  dehnition  in  (IV.49),  we  get 

Dp 


—  +  Po  {dxU  +  dyv)  =  0 

Du  1  „ 

■! - ^xP  =  0 

Dt  Po 

Dv  1  ^ 

H - dyp  =  0 

Dt  Po 

^  +  7Po  {dxU  +  dyv)  =  0  , 


(IV.50) 


(IV.51) 


and  we  can  easily  show  that  this  formnlation  is  eqnivalent  to  a  Lagrangian  scalar 
wave  eqnation 


D^p  _  1P0^2 
DD  po  ^ 

If  we  incorporate  Coriolis  forces,  then  onr  eqnation  set  is 

^  +  Po  {dxU  +  dyv)  =  0 

Du  1  .  „  .  . 

+  —dxP  =  f{v  +  Vo) 
Dt  Po 

—  +  —dyP  =  -f{u  +  Uo) 
JJt  Po 


^  +  7Po  {dxU  +  dyv)  =  0 


(IV.52) 


(IV.53) 
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Using  our  now-familiar  conversion  technique,  we  have 

-  7Po/  {dxV  -  dyu)  (IV.54) 

JJt  Pq 

Now  (IV.49)  and  (IV. 53)  both  imply  ^  so  if  we  apply  ^  to  (IV.54),  we 

get 

^  ~  ^  ~  ^  °  (IV.55) 

Using  the  same  approach  we  used  to  convert  the  no- mean-flow  form  to  the  Klein- 
Gordon  equation,  we  apply  dy  to  (IV. 53b),  to  (IV. 53c),  and  subtract  to  get 

^  {dxV  -  dyu)  =  -f  {dxU  +  dyv)  (IV.56) 


Substituting  this  result  into  (IV.55)  removes  the  vorticity  term  and  gives  us 

_D  /DV  ™  \  Cp  _ 

Dt  \Dt^^  Po  )  Dt 

Integrate  along  characteristics,  assuming  an  initial  value  of  zero,  to  get 

D^p 

V  = 

Po 


Dt^ 


'yP0yj2  f2 
p  = - V  p  —  JP- 


(IV.57) 


(IV.58) 


So  our  system  is  still  equivalent  to  the  Klein-Gordon  equation,  using  Lagrangian 
derivatives  rather  than  time  derivative.  However,  we  will  show  in  Sec.  3  that  such  a 
system  is  inherently  unstable  under  long-time  integration. 


1.  Interior  Discretization  Scheme 

Similar  to  (IV.50),  if  we  dehne  the  Lagrangian  difference 

A  =  At  -f  uqAx  +  voAy  , 
then  the  discrete  basic  system  is  given  by 

Ap  +  Po  {AxU  +  Ayv)  =  0 

V  1  . 

Am  H - AxP  =  0 

Po 

V  1  . 

Av  H - Ayp  =  0 

Po 

Ap  +  7Po  +  Ayv)  =  0  , 


(IV.59) 


(IV.60) 
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which  is  easily  shown  to  be  equivalent  to 

AAp  =  —  {A^A^p  +  AyAyp)  (IV.61) 

Po 

Hence,  the  basic  system  under  advection  should  be  stable  when  implemented  in  con¬ 
junction  with  the  Higdon  NRBC  dehned  for  a  double-sized  grid.  Likewise,  if  we 
include  Coriolis,  we  have  the  discrete  system 


Ap  -f  Po  {A^u  +  Ayv) 

Au  H - AxP 

Po 

Av  +  —AyP 
Po 

Ap  +  7Po  {^xU  +  Ayv) 


0 

/  +  ^  + 

-/  (—p  +  U  +  Uo 
\Po 

0, 


(IV.62) 


and  if  we  use  the  discrete  analog  to  the  continuous  case  derivation  in  the  preceding 
section,  we  can  easily  show  that  this  system  is  equivalent  to 

IPo 


AAp  =  —  {A^A^P  +  Ay  +  Ayp)  -  fp 
Po 


(IV.63) 


which  again  should  be  stable  when  the  Higdon  scheme  is  used  on  a  double-sized 
grid.  Since  the  inclusion  of  gravity  creates  a  system  which  is  not  compatible  in  the 
continuous  case  with  a  wave-like  equation,  its  discrete  formulation  will  also  not  be 
equivalent  to  a  discrete  wave-like  equation.  The  discretization  scheme  in  the  xz  plane 
under  the  influence  of  gravity  is  given  by 

Ap  +  p  {A^u  +  Azw)  =  -p'  {w  -f  Wo) 

Au  -|-  -AxP  =  0 
P 

V  1  A  9 

Aw  + -AzP  =  --p 
P  P 

Ap  +  7p  {A^u  -f  A^w)  =  gp{w^  wo)  ,  (IV.64) 


where  we  replace  the  constants  po  and  po  with  the  ^-dependent  values  p  and  p.  Note 
that  Wo  appears  on  the  right-hand  side  of  the  equations  for  p  and  p  because  of  the 
linearization  process  in  Sec.  H.C.2.  In  practice,  however,  the  presence  of  the  ground 
will  require  Wq  =  0. 
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2.  NRBC  Formulation 

The  addition  of  advection  changes  onr  wave  speed  with  respect  to  the  bonnd- 


ary.  To  determine  the  new  “correct”  wave  speed  estimate,  replace  the  time  derivative 
in  the  Higdon  NRBC  dehnition  (III.  11)  with  the  Lagrangian  derivative  (expressed  for 
a  right-side  bonndary): 


j 

n 


-1-  Cjdx  )  M  =  0 

Dt  ^ 


(IV.65) 


For  waves  striking  the  bonndary,  we  are  only  concerned  with  the  velocity  normal  to  the 
bonndary.  Thns,  we  remove  the  tangential  component  of  the  Lagrangian  derivative, 
expand  and  combine  terms,  and  nse  onr  simplihcation  Cj  =  cq,  giving 


[dt  +  (co  +  Mo)  dxY  u  =  0 


(IV.66) 


Since  onr  interior  discretization  scheme  is  eqnivalent  (or  approximately  so)  to  a  wave 
eqnation  on  a  donble-sized  grid,  onr  NRBC  discretization  for  these  eqnations  is  given 
by 


I  - 
2St 


+  (co  +  Mo) 


2Sx 


j 

a  =  0 


(IV.67) 


3.  Long-Term  Instability  of  Advection  with  Coriolis 

Let  ns  briefly  consider  the  discrete  form  of  (11.63)  in  the  xy  plane  with  Coriolis 
forces.  Using  onr  nsnal  second-order  centered-difference  scheme,  we  have 


Atp  +  uqA^p  +  voAyp  -F  poA^u  +  poAyV 

Atu  +  uqA^u  +  voAyU  +  —A^p 

Po 

Atv  -1-  uqA^v  -F  v^AyV  H - Ayp 

Po 

Atp  +  uqA^p  4-  v^Ayp  +  jpqA^u  +  jpoAyV 


0 

f{v  +  Vo) 

-/(mH-Mo) 

0  (IV.68) 


Assnme  everything  is  initially  at  rest  (no  pertnrbations) ,  and  moMq  ^  0.  With  every 
pertnrbation  variable  eqnal  to  zero  at  times  n  <  0,  onr  valnes  at  time  n  =  1  become 


(IV.69) 


Ki 

2St 


=  fvo 


m  = 

P] 

^  =  0 
2St 

So,  having  started  from  a  zero-perturbation  domain,  the  combination  of  advection 
and  Coriolis  creates  non-zero  perturbations  in  u  and  v.  In  fact,  u  and  v  are  still 
constant  in  the  domain,  but  now  the  constant  is  non-zero.  We  note  that  p  and  p  are 
still  uniformly  zero  throughout  the  domains,  and  it  is  easy  to  show  that  they  will 
always  be  zero.  However,  if  we  look  at  u  and  v  at  the  next  time  step,  noting  that  all 
spatial  differences  are  zero,  we  get 


=  -/K  +  «o 


(IV.  70) 


Substituting  the  previously  determined  values  for  uj,  and  nl  ,  we  get 


ulj  =  2Stfvo-{2Stffuo 

vl  =  -2Stfuo-{2Stffvo  (IV.  71) 

Proceeding  to  time  step  n  =  3,  we  have 

ulj  ^  2  {2Stfvo)  +  O  (^{2Stff^ 

~  -2  (25t/Mo)  +  O  ((25t/)^) 

(IV.  72) 


where  we  are  only  concerned  with  the  terms  which  are  linear  in  2Stf.  It  can  be  shown 
that 


]j  ~  (25^/)  vo  +  0  [{2Stff^ 

j  ^  -[^\{2Stf)uo  +  0{{2Stff)  ,  (IV. 73) 
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where  [-J  denotes  the  floor  function. 

Remark  1:  If  /  7^  0  and  either  uq  ^  0  or  vq  ^  0,  then  our  zero-perturbation 
domain  grows  linearly,  without  bound,  in  u  and/or  v.  While  the  growth  will  be  slow 
due  to  the  magnitude  of  /,  it  will  still  be  present.  Hence,  the  combination  of  Coriolis 
forces  and  non-zero  advection  is  inherently  unstable  in  the  linearized  Euler  equations. 
We  can  only  use  this  equation  set  for  short-time  integrations,  not  for  longer  durations. 

Remark  2:  If  we  perform  a  similar  analysis  in  the  case  of  gravity  and  non¬ 
zero  advection,  we  see  that  that  system  is  inherently  unstable  if  Wq  7^  0.  Horizontal 
advection  does  not  cause  any  problems.  Any  physical  problem  which  considers  gravity 
and  vertical  advection  is  unsuited  to  this  equation  set.  Meteorological  problems, 
however,  will  generally  include  the  ground  and  thus  not  have  any  constant  vertical 
advection. 

4.  Numerical  Examples 

a.  Basic  System,  Infinite  Channel 

For  the  examples  in  this  section,  we  use  the  same  domain  and  initial 
conditions  as  before.  Our  first  example  is  a  basic  system  in  an  infinite  channel  with  a 
constant  lateral  advection  of  Uq  =  100^.  Again  running  the  simulation  up  until  t  = 
24,  we  get  the  error  norms  listed  in  Table  IX.  Fig.  24  shows  the  density  preturbation 
p  for  the  J  =  10  case.  Note  the  faint  wave  propagating  in  the  opposite  direction;  this 
false  wave  is  a  consequence  of  using  centered-differences  for  our  spatial  discretization 
[24].  Fig.  25  compares  the  pressure  perturbation  p  for  low-  and  high-order  NRBCs. 
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J 

E, 

Eu 

Ep 

1 

0.22447 

0.30496 

0.18074 

0.22447 

2 

0.06216 

0.09791 

0.059247 

0.062159 

3 

0.026666 

0.041939 

0.021853 

0.026666 

4 

0.017115 

0.025099 

0.012509 

0.017115 

5 

0.01258 

0.019265 

0.0097557 

0.01258 

6 

0.010007 

0.014986 

0.0074502 

0.010007 

7 

0.0082563 

0.012216 

0.0059657 

0.0082562 

8 

0.0070098 

0.010347 

0.0051175 

0.0070098 

9 

0.0060805 

0.0090438 

0.0044283 

0.0060805 

10 

0.0053682 

0.0079242 

0.003876 

0.0053683 

Table  IX.  Error  norms  with  J  G  1 ...  10  in  an  infinite  channel  with  horizontal  advec- 
tion 


Fignre  24.  Plot  of  p  in  basic  system  (IV.49)  with  left-to-right  advection  with  J  =  10 
in  an  infinite  channel.  (BL)  Compnted  solntion.  (Top)  Reference  solntion;  the  area 
corresponding  to  the  compnted  solntion  is  contained  between  the  vertical  lines.  (BC) 
Reference  solntion  truncated  to  compnted  solntion  domain.  (BR)  Delta  between 
reference  solntion  and  compnted  solntion,  with  error  norm  compnted  by  (IV.  14). 
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Presswe  perturbation,  J*1  Pressure  perturbation.  J=10 


20  40  60  fiO  100  20  40  60  60  100 

Error  norm-0.21817  Error  r>orm-0.0045685 


Figure  25.  Comparison  of  p  in  basic  system  (IV.49)  with  left-to-right  advection 
computed  with  J  =  1  and  J  =  10  in  an  infinite  channel,  with  error  norms  computed 
by  (IV. 14).  (TL)  Computed  solution  for  J  =  1.  (TR)  Computed  solution  for  J  =  10. 
(BL)  Delta  between  reference  solution  and  J  =  1  computed  solution.  (BR)  Delta 
between  reference  solution  and  J  =  10  computed  solution. 
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J 

E, 

Eu 

Ey 

Ep 

1 

0.33816 

0.36566 

0.36566 

0.33819 

2 

0.091059 

0.10563 

0.10563 

0.091066 

3 

0.043514 

0.050827 

0.050827 

0.043517 

4 

0.028633 

0.032601 

0.032601 

0.028635 

5 

0.021446 

0.023794 

0.023794 

0.021447 

6 

0.0172 

0.018581 

0.018581 

0.017201 

7 

0.014495 

0.015508 

0.015508 

0.014496 

8 

0.012603 

0.013347 

0.013347 

0.012604 

9 

0.011967 

0.013011 

0.012734 

0.011969 

10 

0.078775 

0.10216 

0.082933 

0.078757 

Table  X.  Error  norms  with  J  G  1 ...  10  in  an  open  domain  with  diagonal  advection 


b.  Basic  System,  Open  Domain 

Now  consider  the  basic  system  in  an  open  domain  with  a  constant 
diagonal  advection  of  Uq  =  90  m/s,  vq  =  90  m/s.  This  time,  we  get  the  error  norms 
listed  in  Table  X.  Not  snrprisingly,  the  error  norms  for  u  and  v  are  identical  np  to 
J  =  8  (with  V  slightly  better  than  u  in  the  J  =  9  case).  In  this  example,  we  see  that 
the  J  =  10  error  norms  are  larger  than  those  for  J  =  9.  Fig.  26  shows  the  a:-velocity 
pretnrbation  u  for  the  J  =  8  case;  Figs.  27  and  28  shows  the  same  variable  for  the 
J  =  9  and  J  =  10  cases,  respectively.  We  see  from  the  latter  two  that  the  system 
begins  to  destabilize  at  the  top-right  corner.  The  reason  for  this  destabilization  is 
still  nnder  investigation.  Fig.  29  compares  the  |/- velocity  pertnrbation  v  for  low-  and 
high-order  NRBCs. 
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•  0  05 


0 

-0  05 


Figure  26.  Plot  of  u  in  basic  system  (IV. 49)  with  bottom-left-to-top-right  advection 
with  J  =  8  in  an  open  domain.  (TL)  Computed  solution.  (Right)  Reference  solution; 
the  area  corresponding  to  the  computed  solution  is  contained  in  the  center  box.  (CL) 
Reference  solution  truncated  to  computed  solution  domain.  (BL)  Delta  between 
reference  solution  and  computed  solution,  with  error  norm  computed  by  (IV.  14). 


SO  100  -50  0  SO  100  150  200 


Error  t>omv-0  009933$ 


Figure  27.  Plot  of  u  in  basic  system  (IV. 49)  with  bottom-left-to-top-right  advection 
with  J  =  9  in  an  open  domain.  (TL)  Computed  solution.  (Right)  Reference  solution; 
the  area  corresponding  to  the  computed  solution  is  contained  in  the  center  box.  (CL) 
Reference  solution  truncated  to  computed  solution  domain.  (BL)  Delta  between 
reference  solution  and  computed  solution,  with  error  norm  computed  by  (IV.  14). 
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so  100  -50  0  SO  100  150  200 


Error  norm-0  10640 


Figure  28.  Plot  of  u  in  basic  system  (IV. 49)  with  bottom-left-to-top-right  advection 
with  J  =  10  in  an  open  domain.  (TL)  Computed  solution.  (Right)  Reference  solution; 
the  area  corresponding  to  the  computed  solution  is  contained  in  the  center  box.  (CL) 
Reference  solution  truncated  to  computed  solution  domain.  (BL)  Delta  between 
reference  solution  and  computed  solution,  with  error  norm  computed  by  (IV.  14). 


Figure  29.  Comparison  of  v  in  basic  system  (IV.49)  with  bottom-left-to-top-right 
advection  computed  with  J  =  1  and  J  =  8  in  an  open  domain,  with  error  norms 
computed  by  (IV.  14).  (TL)  Computed  solution  for  J  =  1.  (TR)  Computed  solution 
for  J  =  8.  (BL)  Delta  between  reference  solution  and  J  =  1  computed  solution.  (BR) 
Delta  between  reference  solution  and  J  =  8  computed  solution. 


76 


J 

E, 

Eu 

Ey 

Ep 

1 

0.22154 

0.59249 

0.25232 

0.22154 

2 

0.10104 

0.3566 

0.15752 

0.10104 

3 

0.34963 

1.3304 

0.39134 

0.34963 

4 

1.4844 

5.1058 

1.391 

1.4844 

5 

2.3286 

7.4879 

3.3526 

2.3286 

6 

34.219 

98.952 

36.169 

34.219 

7 

71.57 

277.37 

83.65 

71.57 

8 

4971.8 

14265 

4154 

4971.8 

Table  XL  Error  Norms  with  J  G  1 ...  8  in  an  Infinite  Channel  with  Horizontal  Ad- 
vection  and  Coriolis 


c.  Coriolis  Forces,  Infinite  Channel 

Having  shown  already  that  the  system  is  inherently  nnstable  for  long¬ 
time  integrations  involving  both  advection  and  Coriolis,  let  ns  look  briefly  at  the 
effectiveness  of  the  Higdon  NRBCs  for  a  short-dnration  simnlation.  We  first  consider 
an  infinite  channel  with  left-to-right  advection,  Uq  =  100^.  Incorporating  Coriolis 
forces,  we  nse  the  discretization  scheme  (IV. 62)  for  the  interior.  Table  XI  shows  the 
error  norms  for  J  G  1 ...  8.  We  see  some  improvement  between  J  =  1  and  J  =  2,  bnt 
the  system  fails  for  higher  J.  We  will  not  show  any  plots  of  these  resnlts. 
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J 

E, 

Eu 

Ey 

Ep 

1 

0.86499 

0.31392 

0.33843 

0.86506 

2 

0.58828 

0.26882 

0.2833 

0.58833 

3 

6.288 

2.5314 

2.67 

6.2885 

4 

327.16 

134.55 

141.92 

327.18 

5 

14895 

4854.5 

5120.1 

14896 

6 

278880 

87044 

91822 

278900 

Table  XII.  Error  Norms  with  J  G  1 ...  6  in  an  Open  Domain  with  Diagonal  Advection 
and  Coriolis 


d.  Coriolis  Forces,  Open  Domain 

Likewise,  when  we  consider  the  system  in  an  open  domain,  with  the 
NRBC  on  all  fonr  sides,  nsing  a  diagonal  advection,  uq  =  Vq  =  100^,  we  get  the 
resnlts  shown  in  Table  XII.  A  slight  improvement  from  J  =  1  to  J  =  2  is  followed  by 
massive  instability  for  higher  J.  Again,  we  will  omit  any  plots  of  these  resnlts. 
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J 

E, 

Eu 

Ey 

Ep 

1 

0.2251 

0.29999 

0.18206 

0.22442 

2 

0.062733 

0.096051 

0.059583 

0.062517 

3 

0.027088 

0.041293 

0.022049 

0.027048 

4 

0.017524 

0.024613 

0.012558 

0.017488 

5 

0.012913 

0.018886 

0.0098005 

0.012887 

6 

0.010289 

0.014672 

0.0074523 

0.01027 

7 

0.0085063 

0.011944 

0.0059542 

0.0084902 

8 

0.007226 

0.010119 

0.005107 

0.0072123 

9 

0.0062668 

0.0088368 

0.0044166 

0.0062548 

10 

0.0055272 

0.0077406 

0.0038622 

0.0055167 

Table  XIII.  Error  norms  with  J  G  1 ...  10  in  an  infinite  channel  with  horizontal 
advection  and  gravity 


e.  Gravity,  Infinite  Channel 

Tnrning  onr  attention  now  to  the  inclnsion  of  gravity,  we  again  consider 
an  infinite  channel  with  horizontal  advection,  uq  =  100^.  Using  the  discretization 
scheme  (IV. 64),  we  get  the  error  norms  shown  in  Table  XIII.  Unlike  the  Coriolis 
case,  this  time  we  get  good  resnlts  all  the  way  np  to  J  =  10.  Fig.  30  illnstrates  the 
difference  in  pressnre  pertnrbation  resnlts  for  small  and  large  J. 


79 


Pressure  perturbation.  J=1 


20  40  60  80  100 


Pressure  perturbation,  J=10 
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Pressure  perturbation  error.  J=1 


Pressure  pertiabation  error.  J=10 
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Figure  30.  Comparison  of  p  in  gravity-influenced  system  (IV. 64)  with  left-to-right 
advection  computed  with  J  =  1  and  J  =  10  in  an  inhnite  channel,  with  error  norms 
computed  by  (IV.  14).  (TL)  Computed  solution  for  J  =  1.  (TR)  Computed  solution 
for  J  =  10.  (BL)  Delta  between  reference  solution  and  J  =  1  computed  solution. 
(BR)  Delta  between  reference  solution  and  J  =  10  computed  solution. 
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J 

E, 

Eu 

Ep 

1 

0.25917 

0.35427 

0.28085 

0.25639 

2 

0.072395 

0.11036 

0.084279 

0.071655 

3 

0.032031 

0.050496 

0.036455 

0.03179 

4 

0.020847 

0.03178 

0.022262 

0.020674 

5 

0.015383 

0.023689 

0.016319 

0.015259 

6 

0.012198 

0.018435 

0.012387 

0.012104 

7 

0.010073 

0.01515 

0.0099833 

0.0099939 

8 

0.0085547 

0.012836 

0.0084265 

0.0084879 

9 

0.0074187 

0.011229 

0.007244 

0.0073608 

10 

0.0066163 

0.010216 

0.0066934 

0.0065656 

Table  XIV.  Error  norms  with  J  G  1 ...  10  in  a  half-plane  with  horizontal  advection 
and  gravity 


/.  Gravity,  Ground  and  Open  Air 

Looking  now  at  an  “open”  domain  in  the  presence  of  gravity,  we  con¬ 
sider  the  open-air  domain  of  Fig.  22.  This  domain  matches  the  physical  idea  of 
modeling  an  open  atmosphere  bonnded  by  the  gronnd.  While  the  presence  of  the 
gronnd  prevents  ns  from  considering  vertical  or  diagonal  advection,  it  still  presents 
ns  with  a  domain  containing  two  adjacent  open  bonndaries.  We  keep  uq  =  100^  for 
onr  advection  speed  and  rnn  onr  standard  simnlation.  The  error  norms  for  different 
valnes  of  J  are  given  in  Table  XIV.  Fig.  31  illnstrates  the  differences  between  low 
and  high  J. 

F.  SUMMARY 

This  high-order  NRBC  implementation  provides  high  accnracy  with  little  com- 
pntational  overhead.  However,  it  has  three  signihcant  limitations: 

1.  The  formnlation  reqnires  one-sided  high-order  spatial  and  temporal  deriva¬ 
tives. 

2.  These  one-sided  derivatives  limit  the  NRBC  order  to  the  size  of  the  domain. 

3.  Coriolis  and  advection  cannot  be  nsed  together  stably. 
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Figure  31.  Comparison  of  p  in  gravity-influenced  system  (IV. 64)  with  left-to-right 
advection  computed  with  J  =  1  and  J  =  10  in  an  open-air  domain,  with  error  norms 
computed  by  (IV.  14).  (TL)  Computed  solution  for  J  =  1.  (TR)  Computed  solution 
for  J  =  10.  (BL)  Delta  between  reference  solution  and  J  =  1  computed  solution. 
(BR)  Delta  between  reference  solution  and  J  =  10  computed  solution. 


In  the  next  chapter,  we  will  look  at  another  NRBC  implementation,  one  which  does 
not  suffer  from  these  limitations. 
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V.  GIVOLI-NETA  NRBCS  FOR  THE 
LINEARIZED  EULER  EQUATIONS 


A.  INITIAL  IMPLEMENTATION  FOR  FIRST-ORDER 
SYSTEMS 
1.  Derivation 

Let  us  now  consider  the  Givoli-Neta  auxiliary  variable  NRBC  described  in 
Sec.  III.B.2.  That  section  described  this  NRBC  formulation  for  a  scalar  wave-like 
equation.  In  this  chapter  we  apply  it  to  the  linearized  Euler  system.  We  hrst  dehne 
the  auxiliary  variables  in  vector  form: 

(pj  =  (fj+i  (V.l) 

(po  =  (p 


Note  that  we  have  left  off  the  truncation  condition  (pj  =  Q  from  our  formulation. 
As  we  derive  the  implementation  below,  we  shall  see  that  we  need  to  modify  the 
truncation  condition  slightly;  the  modihed  truncation  condition  will  be  given  as  part 
of  the  characteristic-based  derivation  (V.IO). 

Givoli  and  Neta  [39,  40]  show  how  to  manipulate  the  Klein-Gordon  equation  to 
remove  the  normal  derivatives.  We  will  show  briefly  how  to  do  so  with  the  linearized 
Euler  equations.  For  dehniteness,  we  proceed  for  the  right-side  boundary;  thus,  n  =  x. 

Begin  with  the  equation  system  in  vector  form: 

dt0  +  Adx0  +  Bdytp  =  C(p  -f  D  ,  (V.2) 


where 


A  =  ( 

/ 

\ 

T 

kP 

U  V 

P  j 

'  Mo 

Po 

0 

0  ^ 

0 

Uq 

0 

J_ 

A  = 

Po 

0 

0 

Mo 

0 

1  0 

7Po 

0 

Mo  ) 
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^  1^0  0  Po  0  ^ 

0  Vq  0  0 

B  = 

0  0  Uo  - 

^  po 

0  0  7Po  ^'o  / 

In  the  basic  system,  C  and  D  are  both  zero.  With  Coriolis,  we  have 

^  0  0  0  0  ^ 

0  0/0 
0-/00 
^  0  0  0  0  y 

D  =  0  /no  -fuo  0  ) 

With  gravity,  replace  v,  Vq,  Pq,  and  po  in  A  and  B  with  w,  Wq,  p,  and  p,  respectively, 
and  use 


C  = 


D  = 


0  0  -p'  0 

0  0  0  0 
-p4  0  0  0 


0 


0  pp  0  y 


-p'wo  0  0  gpWQ  ) 


Dehne  the  linear  differential  operator  L{(p)  to  be 


L{(p)  =  dt(p  +  Adx(p  +  Bdyp  —  C(p  (V.3) 

By  dehnition,  L{fp)  =  D.  A  more  useful  result  is  the  following  lemma: 

Lemma  V.l  The  auxiliary  variables  <pj  satisfy 

i(ft)  =  0  (V.4) 

for  all  j  G  1 . . .  J.  The  only  exception  to  this  lemma  is  when  both  of  the  following 
conditions  are  met: 
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•  n  =  i/c 

•  ^  0 

Proof.  Apply  the  operator  L  to  (V.l)  to  get 

L((^i)  =  L 

=  dH{L{0))  +  —dt{L{0)) 

=  dfiD  H - dtD 

=  0  (V.5) 


Then  proceed  by  induction  to  prove  the  lemma  for  all  j  □. 

(Note  that  if  both  conditions  in  Lemma  V.l’s  exception  are  met,  then  L  {dz(p)  7^ 
dz  {L  {(p)),  and  we  cannot  proceed  with  this  proof  or  the  overall  derivation.  The  NRBC 
which  handles  these  conditions  will  be  given  in  Sec.  C.) 

Solve  (V.l)  for  and,  using  Lemma  V.l,  substitute  the  result  into  (V.2): 


Aipj+i  =  —A  -  I  dt(pj  -  Bdt(pj  +  CLpj  ,  'ij  el...  J 


(V.6) 


We  now  have  an  auxiliary  variable  equation  which  need  be  dehned  only  on  the  bound¬ 
ary.  However,  as  it  is  dehned  solely  on  the  boundary,  and  given  that  everything  on 
the  boundary  is  initially  zero,  we  have  an  equation  system  which  will  always  be  iden¬ 
tically  zero.  We  will  take  a  characteristic-based  approach  similar  to  that  presented 
in  Sec.  6  of  [57]. 

Let  A  be  the  diagonal  matrix  of  the  eigenvalues  of  A  in  decreasing  order,  with  R 
the  corresponding  right  eigenvectors.  With  two  degrees  of  freedom  for  the  eigenvectors 
associated  with  A  =  mq;  we  choose  eigenvectors  which  do  not  map  one  characteristic 

T 


variable  to  one  state  variable,  which  we  would  get  if  we  chose 

T 


10  0  0 


and 


0  0  10 


for  our  uq  eigenvectors.  Instead,  we  use  the  following: 


A 


Uq 


Uq 


Uq  —  Cq 
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R  = 


(  PO 
CO 

1 
0 

y  poco 
( 0 


1  1 

0  0 

-1  1 


.Po 

CO 


1 

0 


R~^  = 


1 

2 

0 

0 


0  0  -poCo  / 

n  ^  \ 


V 


0^0 


A  =  RAR 


1 

2 

-1 


2poco 
_  1 
2cq 

_  1 


2poco 


(V.7) 


With  this  dehnition,  left-multiply  (V.6)  by  R  ^  and  insert  RR  ^  as  needed  to  get 


A6  =  - /j  +  cel,  +  ^ 

Mi+i  =  (-A-i]adi-Ba,(j  +  c(j ,  j  , 


(V.8) 


where 


0 

B 

C 

D 


R  ^(fj 
R-^BR 
R-^CR 
R-^D 


(V.9) 


However,  we  still  have  a  system  which  is  dehned  exclnsively  on  the  bonndary.  To  over¬ 
come  this  problem,  nse  (V.6)  in  characteristic  form  as  a  natnral  bonndary  condition 
for  the  ontgoing  characteristics.  To  close  the  system,  nse  a  Dirichlet  condition  for  the 
highest-order  anxiliary  variables  corresponding  to  the  incoming  characteristics.  The 
hnal  system  then  is 


iA^ia/o  +  Aei 

Cj  J 


—AdxCo  ~  Bdy^o  +  C^o  +  D 
—Bdy^Q  +  C^o  +  D 
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/--A  =  -BdyC,  +  CCj  ,  j  el...J 


0j  =  0 


(V.IO) 


where  the  subscripts  +  and  —  mean  we  apply  these  equations  only  to  the  outgoing 
and  incoming  characteristics,  respectively;  that  is,  the  hrst  equation  only  applies  to 
the  outgoing  characteristic  variables,  the  second  and  third  equations  apply  to  all  four 
characteristic  auxiliary  variables  (for  each  j),  and  the  hnal  equation  only  applies  to 
the  incoming  characteristic  auxiliary  variables.  For  example,  if  we  have  J  =  1,  Mq  >  0, 
C  =  0,  D  =  0,  and  we  use  standard  hrst-order  forward-differences  in  time,  hrst-order 
backward-differences  in  x,  and  second-order  centered-differences  in  y,  then  (V.IO) 
dehnes  the  matrix  system 


1 

0 

0 

0 

0 

0 

0 

0 

dn 

^a,E,i 

0 

1 

0 

0 

0 

0 

0 

0 

tn 

^b,E,i 

0 

0 

1 

0 

0 

0 

0 

0 

tn 

^c,E,i 

1- A, 

0 

0 

0 

StXa 

0 

0 

0 

^d,E,i 

0 

1  —  Af, 
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0 

0 

StXb 

0 

0 

tn 

^l,a,E,i 

0 

0 

1-Ae 

0 

0 

0 

StXc 

0 

0 
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1  —  Arf 

0 

0 
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StXd 
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^l,c,E,i 

.  0 
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0 

0 

0 

0 
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1  ; 
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\  ^l,d,EA  / 

/  en-i 


n—1 

b,E,i 

cn—1 

'!>c,E,i 


(V.ll) 


87 


where  A  denotes  the  individual  eigenvalues  of  the  diagonal  matrix  A;  the  subscript 
1  denotes  the  hrst  (and  only)  auxiliary  variable;  the  subscripts  a,  b,  c,  d  denote  the 
individual  components  of  the  characteristic/state/auxiliary  variable  vectors;  the  sub¬ 
scripts  E  and  E  —  1  denote  the  eastern  boundary  and  the  point  adjacent  to  it,  respec¬ 
tively;  the  subscripts  i,  i  -|-  1,  and  i  —  1  denote  the  grid  point  in  the  y  direction;  and 
the  superscripts  n  and  n  —  1  denote  the  current  and  previous  time  steps,  respectively. 

Note  that  if  we  are  applying  this  NRBC  in  the  presence  of  gravity,  even  in  a 
horizontal  channel,  we  must  take  care  to  consider  the  ^-derivative  of  the  eigenvector 
matrix  R  when  we  convert  from  state  variables  to  characteristics.  In  such  a  case, 
(V.IO)  becomes 


Ci  / 


dtij  -l-  A^_|_i 

^3  / 


+  {b  R  +  C)  ^0  +  D 

-Bd^io  +  {b  (d,R-^)  R  +  C)i,  +  D 

-ml  +  {b  (d,R-^)  R  +  c)i,  ,  iel...J 


0J  =  o]  _  .  (V.12) 

2.  Incompatibility  with  Zero  Advection 

Note  that,  because  of  the  A  on  the  left-hand  side  of  (V.10b,c),  all  of  the 
eigenvalues  must  be  non-zero.  Requiring  the  advection  be  subsonic  is  trivial;  it  can 
be  a  constraint  of  the  model.  But  the  case  of  zero  advection  is  trickier.  The  answer  is 
to  dehne  a  “false”  non- zero  advection,  small  enough  that  there  will  be  no  drift  from 
one  point  to  an  adjacent  point  over  the  duration  of  the  simulation: 

min(5a:) 


«o  < 


2A 


(V.13) 


where  min(5a:)  is  the  smallest  grid  spacing  in  the  x  direction,  and  tmax  is  the  duration 
of  the  simulation.  If  the  false  Uq  is  too  small,  however,  there  is  a  risk  that  the  resulting 
matrix  could  be  numerically  singular.  Realistically,  this  concern  is  unfounded;  we  have 
seen  experimentally  (using  Matlab)  that  the  matrix  A  is  not  ill-conditioned  so  long 
as  ||mo||  >  (approximately  an  eighth  of  an  inch  per  year). 


J 

Ep 

Eu 

Ey 

Ep 

1 

0.15985 

0.32103 

0.13114 

0.15985 

2 

0.1581 

0.32249 

0.12376 

0.1581 

3 

0.15697 

0.32307 

0.12029 

0.15697 

4 

0.15617 

0.32308 

0.11844 

0.15617 

5 

0.15557 

0.32274 

0.11729 

0.15557 

6 

0.15506 

0.32218 

0.11646 

0.15506 

7 

0.15462 

0.32151 

0.11579 

0.15462 

8 

0.15422 

0.32077 

0.1152 

0.15422 

9 

0.15384 

0.32 

0.11465 

0.15384 

10 

0.15347 

0.3192 

0.11413 

0.15347 

Table  XV.  Error  norms  for  Givoli-Neta  NRBCs  in  basic  system  (IV. 1)  with  J  G 
1 ...  10  in  an  infinite  channel  with  no  advection 


3.  Numerical  Examples 

Let  ns  consider  some  nnmerical  examples.  In  each  example,  we  nse  onr  stan¬ 
dard  domain  and  initial  condition  in  an  infinite  channel  (Fig.  19).  We  hrst  consider 
the  basic  case  with  no  advection.  Since  we  have  to  dehne  a  false  advection  on  the 
bonndary,  we  set 

uo  =  (V.14) 

^^max 

For  this  domain  and  simnlation  dnration,  that  eqnates  to  Uq  =  1.3889^.  Fnrther- 
more,  we  set  Cj  =  1  for  all  j  [23,  55],  a  simplihcation  we  will  nse  for  all  snbseqnent 
derivations  for  these  NRBCs.  Table  XV  shows  the  state  variable  error  norms  for 
J  G  1 ...  10,  Fig.  32  shows  the  state  variable  p  for  the  J  =  10  case,  and  Fig.  33  com¬ 
pares  the  state  variable  u  for  the  J  =  1  and  J  =  10  cases.  We  note  that  there  is  very 
little  difference  between  the  low-order  and  high-order  resnlts.  This  is  a  conseqnence  of 
the  characteristic-based  approach  [72]  as  well  as  a  conseqnence  of  converting  the  nor¬ 
mal  derivative  to  tangential  [36].  If  we  instead  nse  a  horizontal  advection  uq  =  100^, 
we  get  the  resnlts  shown  in  Table  XVI  and  Figs.  34  and  35. 

If  we  incorporate  Coriolis  forces  in  this  domain,  we  see  that  the  Givoli-Neta 
NRBCs  still  prodnce  valid  resnlts.  Tables  XVII  and  XVIII  show  the  error  norms  for 
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Density  perturbation,  reference  solution 
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Figure  32.  Plot  of  p  for  Givoli-Neta  NRBCs  in  basic  system  (IV.  1)  with  J  =  10  in 
an  infinite  channel  with  no  advection.  (BL)  Computed  solution.  (Top)  Reference 
solution;  the  area  corresponding  to  the  computed  solution  is  contained  between  the 
vertical  lines.  (BC)  Reference  solution  truncated  to  computed  solution  domain.  (BR) 
Delta  between  reference  solution  and  computed  solution,  with  error  norm  computed 
by  (IV.  14). 

the  state  variables  for  the  Wq  =  0  and  Uq  =  100  cases,  respectively.  In  the  interest  of 
brevity,  we  will  omit  any  hgures  from  these  examples. 

Similarly,  incorporating  gravity  in  the  xz  plane  poses  no  additional  difficulties 
in  this  horizontal  channel,  as  shown  in  Tables  XIX  and  XX.  We  will  see  in  Sec.  C 
how  to  deal  with  the  effects  of  gravity  on  a  vertical  open  boundary. 

Interestingly,  even  though  these  error  norms  appear  to  decrease  only  slightly 
as  J  increases,  we  can  see  from  a  logarithmic  plot  that  they  do,  in  fact,  decay  ex¬ 
ponentially.  See  Fig.  36  for  a  graphical  representation  of  the  values  in  Table  XX. 
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J 

E, 

Eu 

Ep 

1 

0.1556 

0.23218 

0.15239 

0.1556 

2 

0.15506 

0.23139 

0.15096 

0.15506 

3 

0.15452 

0.2306 

0.14954 

0.15452 

4 

0.15399 

0.22982 

0.14815 

0.15399 

5 

0.15346 

0.22904 

0.14679 

0.15346 

6 

0.15293 

0.22827 

0.14545 

0.15293 

7 

0.1524 

0.22751 

0.14413 

0.1524 

8 

0.15188 

0.22674 

0.14283 

0.15188 

9 

0.15136 

0.22599 

0.14156 

0.15136 

10 

0.15085 

0.22523 

0.1403 

0.15085 

Table  XVI.  Error  norms  for  Givoli-Neta  NRBCs  in  basic  system  (IV. 1)  with  J  G 
1 ...  10  in  an  infinite  channel  with  left-to-right  advection 


J 

E, 

Eu 

Eu 

Ep 

1 

0.12541 

0.25444 

0.1244 

0.12541 

2 

0.12484 

0.25298 

0.11594 

0.12484 

3 

0.12438 

0.25176 

0.11182 

0.12438 

4 

0.12394 

0.25063 

0.10959 

0.12394 

5 

0.12351 

0.24953 

0.10822 

0.12351 

6 

0.12307 

0.24843 

0.10724 

0.12307 

7 

0.12263 

0.24734 

0.10647 

0.12263 

8 

0.12218 

0.24626 

0.10581 

0.12218 

9 

0.12174 

0.24518 

0.10521 

0.12174 

10 

0.1213 

0.24411 

0.10465 

0.1213 

Table  XVII.  Error  norms  for  Givoli-Neta  NRBGs  in  Goriolis-inflnenced  system  (IV.  15) 
with  J  G  1 ...  10  in  an  infinite  channel  with  no  advection 
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J 

E, 

Eu 

Ey 

Ep 

1 

0.24986 

0.6637 

0.26511 

0.24986 

2 

0.24943 

0.66284 

0.26447 

0.24943 

3 

0.24901 

0.662 

0.26385 

0.24901 

4 

0.2486 

0.66119 

0.26325 

0.2486 

5 

0.24819 

0.66039 

0.26268 

0.24819 

6 

0.24779 

0.65961 

0.26213 

0.24779 

7 

0.2474 

0.65886 

0.2616 

0.2474 

8 

0.24702 

0.65812 

0.26109 

0.24702 

9 

0.24664 

0.6574 

0.2606 

0.24664 

10 

0.24626 

0.65669 

0.26013 

0.24627 

Table  XVIII.  Error  norms  for  Givoli-Neta  NRBCs  in  Coriolis-inflnenced  system 
(IV.  15)  with  J  G  1 ...  10  in  an  infinite  channel  with  left-to-right  advection 


J 

E, 

Eu 

Eu 

Ep 

1 

0.12998 

0.25364 

0.1265 

0.1294 

2 

0.1294 

0.25218 

0.11796 

0.12881 

3 

0.12893 

0.25098 

0.11379 

0.12833 

4 

0.12849 

0.24986 

0.11152 

0.12788 

5 

0.12804 

0.24877 

0.11012 

0.12742 

6 

0.12758 

0.24769 

0.10913 

0.12697 

7 

0.12713 

0.24662 

0.10835 

0.12651 

8 

0.12667 

0.24555 

0.10768 

0.12605 

9 

0.12621 

0.24449 

0.10707 

0.1256 

10 

0.12575 

0.24343 

0.1065 

0.12514 

Table  XIX.  Error  norms  for  Givoli-Neta  NRBGs  in  gravity-inflnenced  system  (IV.31) 
with  J  G  1 ...  10  in  an  infinite  channel  with  no  advection 
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X-velocity  perturbation,  J=1 


X-velocity  perturbabon,  J=10 


Figure  33.  Comparison  of  u  for  Givoli-Neta  NRBCs  in  basic  system  (IV. 1)  computed 
with  J  =  1  and  J  =  10  in  an  infinite  channel  with  no  advection,  with  error  norms 
computed  by  (IV.  14).  (TL)  Computed  solution  for  J  =  1.  (TR)  Computed  solution 
for  J  =  10.  (BL)  Delta  between  reference  solution  and  J  =  1  computed  solution. 
(BR)  Delta  between  reference  solution  and  J  =  10  computed  solution. 


J 

E, 

Eu 

Ep 

1 

0.1402 

0.20569 

0.14568 

0.13973 

2 

0.13959 

0.20481 

0.14417 

0.13912 

3 

0.13899 

0.20395 

0.14269 

0.13852 

4 

0.13839 

0.20309 

0.14123 

0.13793 

5 

0.13779 

0.20223 

0.13979 

0.13733 

6 

0.1372 

0.20138 

0.13837 

0.13674 

7 

0.1366 

0.20053 

0.13698 

0.13615 

8 

0.13601 

0.19969 

0.13561 

0.13556 

9 

0.13543 

0.19885 

0.13426 

0.13498 

10 

0.13484 

0.19802 

0.13293 

0.1344 

Table  XX.  Error  norms  for  Givoli-Neta  NRBCs  in  gravity-influenced  system  (IV.31) 
with  J  G  1 ...  10  in  an  infinite  channel  with  left-to-right  advection 
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X-direction  velocity  perturtation.  reference  solution 


•50  0  50  100  150  200 

X-direction  velocity  perturbation  Truncated  reference  solution  Solution  delta 


50  100  SO  100  50  100 

Error  norm-0.19874 


Figure  34.  Plot  of  u  for  Givoli-Neta  NRBCs  in  basic  system  (IV. 1)  with  J  =  10 
in  an  infinite  channel  with  left-to-right  advection.  (BL)  Computed  solution.  (Top) 
Reference  solution;  the  area  corresponding  to  the  computed  solution  is  contained 
between  the  vertical  lines.  (BC)  Reference  solution  truncated  to  computed  solution 
domain.  (BR)  Delta  between  reference  solution  and  computed  solution,  with  error 
norm  computed  by  (IV.  14). 
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Figure  35.  Comparison  of  p  for  Givoli-Neta  NRBCs  in  basic  system  (IV.  1)  computed 
with  J  =  1  and  J  =  10  in  an  infinite  channel  with  left-to-right  advection,  with  error 
norms  computed  by  (IV.  14).  (TL)  Computed  solution  for  J  =  1.  (TR)  Computed 
solution  for  J  =  10.  (BL)  Delta  between  reference  solution  and  J  =  1  computed 
solution.  (BR)  Delta  between  reference  solution  and  J  =  10  computed  solution. 


Density  perturbation  error  norm 


X-velocity  perturbation  error  norm 


Z-velocity  perturbation  error  norm  Pressure  perturbation  error  norm 


Figure  36.  Logarithmic  plot  of  state  variable  error  norms  (IV.  14)  for  Givoli-Neta 
NRBCs  applied  to  gravity-influenced  system  (IV. 31)  with  J  G  1 ...  10  in  an  inhnite 
channel  with  left-to-right  advection.  (TL)  Error  norms  for  p.  (TR)  Error  norms  for 
u.  (BL)  Error  norms  for  w.  (BR)  Error  norms  for  p. 
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B.  CORNER  CONDITIONS  IN  AN  OPEN  DOMAIN 
1.  Derivation 

In  an  anxiliary  variable  NRBC  scheme,  the  corner  where  two  open  bonndaries 
intersect  is  a  sonrce  of  concern.  Assnming  the  tangential  derivatives  are  approximated 
with  a  centered-difference  formula,  points  adjacent  to  the  corner  depend  on  the  value 
at  the  corner.  (This  is  not  the  case  in  a  standard  high-order  NRBC  such  as  the 
Higdon  sequence,  where  there  are  no  tangential  derivatives  and  thus  no  dependence 
on  the  corner  values,  as  we  noted  in  Sec.  IV. B. 4.)  So  how  should  we  treat  this  corner? 
Is  it  to  be  considered  part  of  one  boundary  but  not  the  other?  Even  so,  we  still  have 
to  use  a  different  discretization  scheme  than  that  used  elsewhere  on  the  boundary. 
We  propose  here  to  treat  the  corner  as  if  it  were  part  of  a  curved  boundary;  dehne  the 
normal  derivative  as  the  45°  vector  bisecting  the  normal  vector  of  the  two  adjacent 
sides  (see  Fig.  37).  With  this  dehnition,  we  have  for  the  top-left  and  top-right  corners 


Figure  37.  The  “normal”  derivative  for  the  top- right  corner  of  an  open  domain 


of  Fig.  22, 


Top-left:  dn 
Top- right:  dn 


-dx  +  d^ 

V2 

dx  +  d^ 

V2  ’ 


(V.15) 
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and  the  associated  tangential  derivatives  are  given  by 


Top- left:  dr 
Top- right:  dr 


dx  +  d^ 

72 

dx  -  d^ 

72 


(V.16) 


In  order  to  derive  onr  characteristic  system  for  the  anxiliary  variables,  we  mnst  rewrite 
(V.2)  in  terms  of  normal  and  tangential  derivatives.  It  is  straightforward  to  show  that 


dt(p  +  Ndn0  +  Tdr0  =  C(p  -f  D  , 


(V.17) 


where 


Top-left:  N 
T 

Top- right:  N 
T 


1 

72 

1 

72 

1 

72 

1 


i-A  +  B) 
(A  +  B) 
(A  +  B) 
(A-B)  . 


(V.18) 


With  these  dehnitions,  it  is  easy  to  show  that  the  characteristic  NRBC  system  for 
each  corner  is  given  by 


'dtCo  =  -Ad,Co-fdrCo  +  CCo  +  D]^ 

(/-A)at^+Ag+1  =  -fdrCj  +  CC,  (V.19) 

[6  =  o]_, 

where 

A  =  diag  ^  no  -f  Co  no  no  no  -  Cq  ^ 

A  =  RNR-^ 

f  =  R-^TR  ,  (V.20) 


and  C  and  D  dehned  as  in  (V.9).  On  the  top-left  corner,  we  have 


To  = 


Uq  +  Wq 

72 


R  = 


R-^  = 


CQ 


\/2po 


1  -1 


\/2po 

_CQ_ 


V 


0  0 

\/2po 


CQ 

\/2po 

CQ 

\/2po 


4co 

1 

4 

1 

4 

\/2po 

4co 


\/2po 

4co 

1 

4 

1 

4 

_  \/2po 
4co 


2c; 


2c; 


and  on  the  top-right  corner,  we  nse 


no 

t-q 


Wo  +  Wo 

72 

Wo  -  Wo 


72 


/ 


CQ 

\/2po 

CQ 

\/2po 


''0 


-1 


n  \/2po 
4co 


1 

4 

_  1 
4 

\/2po 

4co 


1  1  \ 

_1  CQ 

\/2po 

1  CQ 

\/2po 

0  eg  J 

V^PO  1  \ 
4co  2c^  I 


1  1 


\/2po _ 1  I 

4co  ^  / 


In  implementing  the  hnite-difference  approximations  of  5^-,  decompose  it  into  its  con- 
stitnent  dx  and  dz  components.  Use  the  hr  anxiliary  variable  for  dx,  and  nse  the  F l 
or  Fr  anxiliary  variable  for  dz,  nsing  one-sided  differences.  Note  also  that  the  eigen- 
valnes  of  N  mnst  be  non-zero.  We  can  again  nse  a  “false”  advection  to  get  aronnd 
this  difficnlty,  if  needed. 
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J 

E, 

Eu 

Ey 

Ep 

1 

0.16986 

0.20783 

0.19382 

0.16986 

2 

0.16904 

0.20686 

0.19168 

0.16905 

3 

0.16823 

0.20591 

0.18957 

0.16823 

4 

0.16742 

0.20496 

0.18749 

0.16742 

5 

0.16661 

0.20401 

0.18545 

0.16661 

6 

0.16581 

0.20307 

0.18344 

0.16581 

7 

0.16501 

0.20214 

0.18146 

0.16502 

8 

0.16422 

0.20121 

0.17951 

0.16422 

9 

0.16343 

0.20029 

0.1776 

0.16343 

10 

0.16265 

0.19938 

0.17571 

0.16265 

Table  XXL  Error  norms  for  Givoli-Neta  NRBCs  in  basic  system  (IV. 1)  with  J  G 
1 ...  10  in  an  open  half-plane  with  left-to-right  advection 

2.  Numerical  Examples 

For  simplicity,  let  ns  consider  jnst  one  example  of  this  open-air  domain.  We 
nse  onr  standard  example  with  left-to-right  advection  with  no  inhomogeneons  forces. 
We  reqnire  a  “false”  vertical  advection  on  the  top  bonndary  and  at  the  corners;  we 
nse  Wq  =  =  1.0417^.  Table  XXI  and  Fig.  38  show  the  error  norms  for  the 

state  variables  for  J  G  1 ...  10.  Again,  althongh  the  differences  are  small,  they  are 
nonetheless  exponential. 

C.  GRAVITATIONAL  EFFECTS 
1.  Derivation 

Onr  claim  that  L{(pj)  =  0  in  the  preceding  section  hinges  on  the  fact  that  the 
coefficient  matrices  are  constant  in  the  direction  of  the  normal  derivative.  However, 
this  is  not  the  case  when  we  consider  an  NRBC  on  Ft  in  the  presence  of  gravity.  In 
that  case,  the  coefficient  terms  A,  B,  C,  and  D  all  depend  on  Recall  the  differential 
operator  L[(p)  dehned  in  (V.3).  We  still  have  L{(p)  =  D  \yy  dehnition,  bnt  what  of 
the  anxiliary  variables?  Since  L  {dz(p)  7^  {L  {(p)),  we  come  np  against  the  exception 
to  Lemma  V.l.  A  new  derivation  is  needed.  Apply  the  operator  L  to  both  sides  of 
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Figure  38.  Logarithmic  plot  of  state  variable  error  norms  (IV.  14)  for  Givoli-Neta 
NRBCs  applied  to  basic  system  (IV. 1)  with  J  G  1 ...  10  in  an  open  half-plane  with 
left-to-right  advection.  (TL)  Error  norms  for  p.  (TR)  Error  norms  for  u.  (BL)  Error 
norms  for  v.  (BR)  Error  norms  for  p. 


(V.la,  i  =  0),  and  simplify 


i(^i) 


i(^i) 


L(d,(^  +  L(d,0) 

/ 

\dttp  +  Ad^t<^  -  Cdt(p\ 

< 

+  [dzS  +  Ad^^(p  -  Cd^(p] 

/ 

dt  (dtp  +  Ad^p  +  Bd^p  -  Cp) 

<  +d^  (dtp  +  Ad^p  +  Bd^p  -  Cp) 

-(d^A)d^p-  (d^B)d^p  +  (d^C)p 
dt(L(p))  +  d^(L(p))  -  (d^A)d^p-  (d^B)d^p  +  (d^C)p 

dtD  +  d^D  -  (d^A)d^p-  (d^B)dzp  +  (d^C)p 

-(d^A)d^p-  (d^B)d^p  +  (d^C)p  +  d^D  (V.21) 


As  we  continue  applying  the  operator  L  to  successive  Pj's,  this  expression  will  result 
in  ever-higher  2;-derivative  terms,  exactly  the  case  we  are  trying  to  avoid.  We  need 
to  dehne  our  reference  states  in  such  a  way  as  to  remove  this  difficulty.  The  most 
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effective  approach  is  to  use  the  exponentially-decaying  atmospheric  model  (IV. 32), 
given  in  Sec.  IV.D.l. 

We  can  rewrite  the  matrices  A  and  B  as  the  sum  of  a  constant  matrix  and  a 
^-dependent  one: 


A  —  rto/  +  Az 

B  =  wqI  Bz  ,  (V.22) 

where  Az  and  Bz  are  the  ^-dependent  cells  of  A  and  B,  respectively.  Given  these 
^-dependencies,  and  the  dehnitions  for  p  and  p,  we  have 


dzA  =  —aAz 
dzB  =  —aBz 
dzC  =  -aC 
dzD  =  -aD  . 


Continuing  our  derivation  of  L{(pi)  using  these  dehnitions,  we  have 


L{ipi)  =  aAzdx0  +  aBzdzff  —  aC(p  —  aD 

{l{(p)  -  D  -  dt(p  -  UQdx(p  -  wodz^ 
{p  -  D  -  dtif  -  uodx<f  -  wodz(p 
=  -a  {dt0  +  UQdx0  +  wodz(p) 

L{0i)  =  -adF0 , 


=  a 


=  a 


(V.23) 


(V.24) 


where  dp  denotes  the  Lagrangian  derivative  along  the  how.  If  we  continue  this  deriva¬ 
tion  for  successive  (pj  terms,  we  get 

L{02)  =  —2adFPi  —  a^dpp 

Lips)  =  —3adFP2  —  Sai^dFPi  —  a^dpp  (V.25) 

etc. 


This  leads  us  to  the  following  lemma: 
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Lemma  V.2 


k=l 


(V.26) 


Proof.  Having  demonstrated  it  for  the  j  =  1  case,  we  now  proceed  by  indnction: 


=  < 


L{(pj+i)  =  L{dt(pj)  +  L{d^(pj) 

Ozt^j  “I”  -^dxz^j  “1“  BOzz^j  COz^j 
dt{L{(fj))  +  dz  +  Adx0j  +  Bdz<fj  -  CCpj) 

-{dzA)dx<^j  -  {dzB)dz(fij  +  {dzC)(pj 
dt{L{0j))  +  dz{L{0j))  +  aAzdx0j  +  aBzdz0j  -  aCtpj 
dt{L{^j))  +  dz{L{^j)) 

+a(L((pj)  -  dt(pj  -  uodx^Pj  -  wodz<fj) 

+  dz{L{(fj))  +  aL{(pj)  -  adp^j 

dt  (-  ELi  {i)a’^dp<fj-k)  +  dz  (-  Ei=i  {i)a’^dp<fj-k) 
{  +a  (-  Ei=i  {l)a'^dp(fj_k)  -  adp(fj 

.fco  f a  I  o  ,7?  ^  (A  ^k+l; 


=  < 


=  < 


(V.27a) 

(V.27b) 

(V.27c) 

(V.27d) 

(V.27e) 

(V.27f) 

(V.27g) 


“  Yj  {dt0j-k  +  dz0j-k)  -  Y  [A<y^""^dp0j-k  (v.27h) 

fc— 1  \  /  /c^O  \  / 


j 

-i: 

k=l 

^  O'  +  1 


i+1 


dpLpj+i-k  -Y^{^j^^_ij  Oi^dp(fj+i-k 
a^dp(pj+i_k  -  a^^'^dpip 


2  J 

k  [k-l 


Y 

k=l 


k 


a 


'dp(pj+i_k  -  a^^^dp0 


E  Pt  ° 

k=i  '  ^ 


(V.27i) 

(V.27j) 

(V.27k) 

(V.28) 


We  can  now  nse  this  L{(p)  to  remove  the  normal  derivative  terms.  Proceeding  with 
each  (pj  in  snccession,  we  hrst  have 


Bdz<p  =  —dt0  —  Adx(p  +  C(f  +  D  . 


(V.29) 


Next,  we  have 


L{(pi)  =  dt0i  +  Adx(fi  +  Bdzifi  -  C(fi  =  -adp(fi  , 


(V.30) 
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which  we  combine  with  (V.29)  to  get 


Bdz0i  =  -  Adx(fi  +  C0i  -  adt(p  -  au^d^Cp  -  awodz(p 

/ 

-dt(pi  -  Adx(pi  +  C(pi  -  a  {I  -  woB~^)  dt(p 
—a{uQl  —  WQB~^A)dx(p  —  oiWQB~^{C(p  +  D^  . 


Lemma  V.3 


Bd^ifj 


-dt(pj  -  Ad^ifj  +  Cipj  -  Ei=i  {Cjoi'^Pk  [(/  -  wqB  dt(pj-k 
<  +  {uqI  -  wqB-^A)  d^0j_k  +  woB~^C0j_k] 

—a^ PjWoB~^D  , 


for  all  j  >  1,  where  the  polynomial  sequence  Pj  is  defined  recursively  by 


Pj  {woB-^)  =  ^  -  E  {woB-^)  . 

Proof.  As  with  the  preceding  lemma,  having  demonstrated  the  j  =  1 
proceed  by  induction: 


+  Adx^fj+i  +  Bdz^Pj+i  —  Cifjj^i 


—  Adx^j^i  +  C(pjj^i  —  ^ 


i+l  +  IN 


(y^dpfij+i-k 


^’(A+i) 
t+i 

-  E 

k=l 


k=l 


i  1  _  fc 


k 


a  {dt  +  uodx  +  wqOz)  fij+i-k 


<h((Pj+i)  -  ELi  {dtfij+i-k  +  uodxfij+i-k 

PwqB~^  {—dtfij+i-k  —  Adxfij+i-k  +  Cfij+i-k 
-  eS’^'  [(/  -  woB~^)  dtfi^+^_k-i 

+  {uqI  -  wqB-^A)  dxfij+i-k-i  +  woB-^Cfij+i_k-i] 

-a^+^-^Pj+i_kWoB^^D^'^  -  {dtfi  +  uodxfi 
+woB~^  (^-dtfi -  Adxfi  +  Cfi  + 


(V.31) 

(V.32) 

(V.33) 
case,  we 

(V.34) 

(V.35a) 

(V.35b) 

(V.35c) 
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Bd^iPj+i 


'  »(»5i+i)  -  Ei=i  [(/  -  »oS^')  d,^Pj+i-k 

+  {uol  -  woB~^A)  da,(pj+i-k  +  woB-^C(fj+i_k\ 

+»oB^‘  Ei.i  ec;-'‘  A]-'‘ypi  k/  - 

+  {uqI  -  wqB^^A)  dx0j+i-k-i  +  woB~^C0j+i-k-i] 

+woB~^  ELi  {^l^)a^^^Pj+,.kWoB-^D  -  ((/  -  woB'^)  dt^p 

+  {uqI  —  wqB~^A)  dxP  +  wqB~^  {C(p  + 

^  $(^,+i)  -  ES  [(/  -  WoB-^)  dt0,+,.u 

+  {uqI  -  wqB-^A)  d^0j+i_k  +  woB^'^Cpj+i.k] 

<  +WoB-^  ELi  ELL'  (L')  [(^  -  ^oB-^)  dtL+ilY^?^) 

+  (uqI  -  woB~^A)  dxPj+i-k-i  +  WQB~^C(pj+i-k-i] 


[  -a^+i  (/  -  ELi  WoB-^D  . 

Since  the  binomial  coefficients  are  symmetric,  we  can  replace  Pj+i-k  with  Pk  in  the 
last  line  of  the  previons  eqnation.  The  terms  in  the  hnal  gronp  of  parentheses  are 
thns  Bj+i  by  dehnition.  As  a  temporary  shorthand,  let 


^Lj+l-k-l) 

n 


(/  -  wqB  dtPj+i-k-i 

< 

+  {uqI  -  wqB^^A)  d^pj+i-k-i  +  woB-^Opj+i-k-i 

L^^Pj+iw^B'^D  . 


This  gives  ns 


Bdz^j-\-i  \ 


<i>(^i+i)  -  Eis  (':>'•  [(/  -  woB-')  m+i-, 

+  {uqI  -  woB~^A)  d^(pj+i_k  +  woB~^C(pj+i_k] 

+i»oB->  Ei.,  E&y‘  (y‘)  ('+i-‘)a‘*'-P'®(^j+i-‘-o 


(V.37) 


n 


We  wonld  like  to  replace  the  donble  snmmation  with  one  in  terms  of  rather  than 
Let  r  =  k  +  1.  Looking  at  the  valnes  assnmed  hy  k  +  I  and  /,  and  recognizing 
that  I  <  k  +  I  hy  dehnition,  we  can  rewrite  this  donble  snmmation  as 


r=l l=l 


I 
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A  quick  check  shows  that,  not  only  are  the  indices  mapped  correctly,  we  also  have 
the  correct  number  of  terms  in  the  new  summation: 


j  j+l-k  j+lr-l  .V.-il'i 

E  E  i  =  EEi  =  '-^ 

k=l  1=1  r=l 1=1  ^ 


Also,  we  can  simplify  the  binomial  expansions  by 


O'  + 


Returning  to  our  derivation,  and  reverting  r  back  to  k,  we  have 


Bd^(pj+i  =  I 


=  < 


=  < 


-  ES  !(/  - 

+  {uqI  -  WqB-^A)  +  WoB-^Cipj+i-k] 

[  eS  Eti‘  (d‘)  (?)a‘-Pie(ft+i-<=)  -  » 
Hn+i)  -  ntt 

[  eS  Eti‘  (d‘)  (?)«‘-Pie(ft+i-t)  -  Si 
Eili  (/  -  Eti‘  0)»>oB-‘P,) 


V  k 
i+1 


J  +  1  \  ^  fc  ; 


k=l 


k 


Bd^ifj+i  = 


^tAj+i  A  C(pjj^i 

-  e£i  ^x^y^Pk  [(/  -  woB~^) 

+  {uqI  -  wqB^^A)  dx0j+i-k  +  woB~^C0j+i-k] 


-a- 


’+^Pj+iWoB-^D  □ 


(V.38a) 

(V.38b) 

(V.38c) 

n 

(V.38d) 

(V.39) 


We  are  now  able  to  remove  the  normal  derivative  terms  from  the  auxiliary  variable 
equations.  The  auxiliary  variable  formulation  for  Fy  is 


0j+i  =  {dt  +  0j  (V.40) 

<fo  =  <f 

0j  =  0  . 
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Using  our  formula  for  dz(pj,  we  have 


dt(pj  +  B  ^  {-dt(pj  -  Ada^ipj  +  C0j 

-  ELi  (i)a*U  l(^  -  k'oB'’) 

+  (uqI  -  woB~^A)  dx(pj-k  +  WQB~^C(pj-k] 
—a^ Pj'WQB~^D^  . 


(V.41) 


Collecting  all  dt  terms  on  the  left-hand  side  with  the  (pj+i  term,  leaving  all  the  other 

terms  on  the  right-hand  side,  and  left-multiplying  everything  by  5,  we  have 

/ 

/  ■\  — Adx^-Pj  d"  C(pj 

Y:l-i(i)a'"PkU-woB~Adt0.k  -A. 

^  =  -j:i=AiyPk[{uoI-woB^^A)dx0,.k 

+  {I-B)  dt0i  +  B0.+1  ^ 

+woB~^C0j-k]  -  a^PjWoB~W  . 

(V.42) 

In  characteristic  form,  we  left-multiply  both  sides  by  R~^  and  insert  as  needed 


to  get 

Ei.i 

+  (/  —  A)  dtij  +  A^_|_i 

for  j  G  1,  2, . . . ,  J  —  1,  where 


—Adxij  +  Cij 

-  ELi  [{uoi  -  woA-^a)  dxij-k 

+WQK~^Cij-k  -  a^PjWoA~^D  , 

(V.43) 
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A  =  R-^AR 

P,  =  P-^P,P  (V.44) 

and  ^j,  C,  and  D  are  defined  as  in  (V.9).  As  with  the  simple  case,  we  nse  natnral 
bonndary  conditions  for  the  ontgoing  characteristics  and  the  trnncation  condition  for 
the  incoming.  Hence,  onr  anxiliary  variable  system  is 


(/-A)a,6  +  Aei 


-Ad^io  -  Ad^io  +  {A{d,R-^)R 
~Adx^o  +  +  i) 


E 


3 

k=l 


\  -Adxl  +  Cl 


g)aVa/-woA-i)aV_fc  1 

>  =  < 

-EJ.1 

+  (/  —  A)  dtij  +  A^-+i  j 

-WqA-^  A)dxij-k 

^  AwqA 

- 

[e 

/  =  0 

1 

5 

{uqI 


a^PjWoA 

(V.45) 


where  we  again  nse  the  +  and  —  snbscripts  to  denote  ontgoing  and  incoming  char¬ 
acteristics,  respectively.  Note  also  that,  in  order  to  have  the  term  in  the  first 
eqnation,  we  nse  the  prodnct  rnle  and  the  ^-dependence  of  R~^: 

dz(o  =  dz  (P^  Vo)  =  R~^d^(pQ  +  {d^R~^^  (po  . 

In  the  absence  of  gravity,  a  =  0,  and  (V.45)  is  analogons  to  (V.IO). 

We  lose  some  compntational  efficiency  with  this  implementation.  Eq.  (V.IO) 
defined  a  bi-diagonal  matrix  system,  which  reqnires  0{J)  operations  per  time  step 
with  an  efficient  matrix  solver.  On  the  other  hand,  (V.45)  is  lower-Hessenberg  and 
reqnires  0{J^)  operations  per  time  step. 

When  we  extend  this  derivation  to  the  open  domains  defined  in  the  previons 
section,  we  get  the  following  system  nsing  N  and  T  instead  of  A  and  B  (compare  to 
(V.19)): 

-AaV-  fdxi  +  (A  {dnR-^)  R 

+f{drR-^)R  +  C)iAD\^ 
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{I-A)dtio  +  Hi  =  -f[dM-{drR~^)Rio)+C(o  +  D 

[  -f{d^(j-{d^R-^)RQ+C(j 


ELi  (/  -  noA-i)  dtCj-k 
+  (/  —  A)  dtCj  +  A0j+i 


>  =  < 


0 


=  0 


)  (Brij-t  -  (drR-')  Rij-k] 

[  +noA~^C^j-k  -  a^PjnoK~^b 

(V.46) 


where 


Pk  =  Pk{noN-^) 


A  —  diag  no  +  Co  no  uq  no  -  co 

N  =  R-^AR 


f  =  R-^TR 


a 


a  = 


V2  • 


In  addition,  we  must  replace  po  in  the  eigenvector  matrices  R  and  R~^  with  the 
dependent  p. 

About  the  Polynomial  Sequence.  A  quick  literature  search  [2,  60,  103,  118] 
revealed  no  matches  to  the  polynomial  sequence 


n—1 


Pn{x)  =  1  -  Xj  • 


The  hrst  hve  polynomials  of  this  sequence  are 


Pi  = 

1 

P2  = 

1  -2a: 

Ps  = 

1  —  6a:  +  6a:^ 

Pa  = 

1  —  14a:  +  36a:^  - 

-24a:^ 

P5  = 

1  -  30a:  +  ISOa:^ 

-  240a:2  +  120a:^ 

The  hrst  three  are  constant  multiples  of  the  Bernoulli  polynomials  [103],  but  not  the 
subsequent  ones.  We  have  not  found  any  other  matches  to  this  sequence.  However, 


108 


Dr.  Pante  Stanica  at  the  Naval  Postgraduate  School  observed  [104]  that  the  coeffi¬ 
cients  of  Pn  can  be  generated  directly,  without  use  of  the  recursion  sequence,  by  the 
formula 

where  |^|  denotes  the  Stirling  number  of  the  second  kind  [46].  Other  properties  of 
these  polynomials,  such  as  orthogonality  and  a  generating  function,  may  be  explored 
in  future  research. 

2.  Numerical  Examples 

Let  us  consider  a  few  numerical  examples,  using  our  standard  simulation. 
First,  we  look  briefly  at  the  system  in  a  semi-inhnite  channel,  open  on  top  (see 
Fig.  6),  subject  to  gravitational  forces.  There  is  no  real  advection,  but  we  dehne  our 
false  boundary  advection  wq  =  =  1.0417^  to  keep  the  matrices  non-singular. 

Table  XXII  shows  the  error  norms  resulting  from  the  cases  J  G  1 ...  10. 


J 

E, 

Eu 

Ev 

Ep 

1 

0.073179 

0.08778 

0.11305 

0.049617 

2 

0.073086 

0.087435 

0.11249 

0.04944 

3 

0.07298 

0.08704 

0.11197 

0.049261 

4 

0.072874 

0.086646 

0.11146 

0.049083 

5 

0.07277 

0.086255 

0.11095 

0.048907 

6 

0.072667 

0.085866 

0.11045 

0.048732 

7 

0.072565 

0.08548 

0.10995 

0.048558 

8 

0.072464 

0.085096 

0.10945 

0.048386 

9 

0.072364 

0.084716 

0.10896 

0.048215 

10 

0.072264 

0.084337 

0.10848 

0.048046 

Som. 

0.092721 

0.10278 

0.22495 

0.092238 

Table  XXII.  Error  norms  (IV.  14)  for  Givoli-Neta  NRBCs,  J  G  1 ...  10,  gravity- 
influenced  system  (IV. 31),  semi-inhnite  vertical  channel,  no  advection.  The  error 
norm  from  using  a  Sommerfeld  boundary  condition  is  included  for  comparison. 


Next,  we  consider  the  open-air  domain  (Fig.  22).  We  begin  hrst  with  “zero” 
advection  (false  boundary  advection  uq  =  =  1.3889^,  wq  =  =  1.0417^, 
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chosen  so  that  uq  7^  Wq  to  keep  N  non-singular).  Fig.  39  shows  the  state  variable  p 
at  the  end  of  the  run  for  J  =  10.  Table  XXIII  shows  the  error  norms  (IV.  14)  for  each 


Pressure  perturtation,  reference  solution 


Pressure  perturbation 
100 
60 
60 
40 
20 


Tnjncated  reference  solution 

100  I - 
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ISO  200 

Solution  deta 


i 


Figure  39.  The  solution  for  the  pressure  perturbation  p  using  J  =  10,  Givoli-Neta 
NRBCs  in  an  open  half-plane  subject  to  gravitational  forces  with  no  advection. 
(Top)  Reference  solution;  the  area  corresponding  to  the  computed  solution  is  con¬ 
tained  in  the  bottom-center  box.  (BL)  Computed  solution  using  NRBCs.  (BC)  Ref¬ 
erence  solution  domain  corresponding  to  NRBC  solution  domain.  (BR)  Delta  between 
computed  and  reference  solutions,  with  error  norm  computed  by  (IV.  14). 


state  variable  as  J  goes  from  1  to  10. 

For  our  third  example,  we  use  the  non-zero  mean  flow  Uq  =  100^  (maintaining 
the  same  “false”  Wq  on  the  boundary),  we  get  the  error  norms  in  Table  XXIV. 

In  both  cases,  we  observe  that  there  is  very  little  improvement  when  we  in¬ 
crease  J.  This  property  of  auxiliary  variable  NRBCs  has  been  noted  before  in  other 
contexts  [36]  and  is  the  result  of  discretization  errors  induced  by  converting  the  nor¬ 
mal  derivative  terms  to  tangential.  Higher-order  discretization  schemes  may  improve 
these  results  [55],  but  for  the  purpose  of  this  dissertation,  it  is  sufficient  to  demon¬ 
strate  how  to  use  the  auxiliary  variable  NRBC  methods  with  the  linearized  Euler 
equations.  Furthermore,  we  note  by  the  bottom  row  of  Tables  XXII-XXIV  that  the 
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J 

E, 

Eu 

Ey 

Ep 

1 

0.30742 

0.48811 

0.20486 

0.21176 

2 

0.30615 

0.48976 

0.20039 

0.20991 

3 

0.30532 

0.49034 

0.19818 

0.2087 

4 

0.30472 

0.49024 

0.19688 

0.20784 

5 

0.30426 

0.48972 

0.19598 

0.20717 

6 

0.30387 

0.48896 

0.19524 

0.20662 

7 

0.30353 

0.48805 

0.19459 

0.20613 

8 

0.30321 

0.48706 

0.19399 

0.20568 

9 

0.30291 

0.48602 

0.19342 

0.20526 

10 

0.30261 

0.48497 

0.19286 

0.20485 

Som. 

0.30224 

0.85386 

0.30744 

0.29872 

Table  XXIII.  Error  norms  (IV.  14)  for  Givoli-Neta  NRBCs,  J  G  1 ...  10,  gravity- 
influenced  system  (IV.31),  open-air  domain,  no  advection.  The  error  norm  from 
using  a  Sommerfeld  boundary  condition  is  included  for  comparison. 

auxiliary  variable  method  is  signihcantly  more  accurate  than  the  basic  Sommerfeld 
(hrst-order  Higdon)  condition. 

Note:  Examination  of  our  characteristic  systems  (V.45)  and  (V.19)  reveals 
that,  in  a  uniform  domain  (zero  perturbations  in  all  state  variables),  the  presence 
of  the  non-zero  D  on  the  right-hand  side  will  generate  non-zero  results  at  the  next 
time  step.  Physically,  this  should  not  happen;  it  is  a  consequence  of  replacing  the 
normal  derivatives  with  tangential  derivatives.  To  avoid  this  problem,  a  check  should 
be  made  at  each  point  on  the  boundary  prior  to  applying  the  boundary  condition:  If 
the  solution  at  the  point,  its  immediate  neighbors,  and  the  corresponding  auxiliary 
variables  are  all  zero,  then  the  solution  at  the  point  should  remain  zero. 

D.  SUMMARY 

In  this  chapter  we  have  modihed  the  Givoli-Neta  auxiliary  variable  method  for 
implementation  in  a  hrst-order  system.  Although  our  original  context  is  the  linearized 
Euler  equations,  the  linear  matrix  equation  can  be  extended  easily  to  any  hrst-order 
system.  The  gravity-inhuenced  derivation  is  specihc  to  the  linearized  Euler  equations. 
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J 

E, 

Eu 

Ey 

Ep 

1 

0.20976 

0.25767 

0.24514 

0.18111 

2 

0.20923 

0.25683 

0.24286 

0.18049 

3 

0.20866 

0.256 

0.24061 

0.17987 

4 

0.20813 

0.25517 

0.2384 

0.17925 

5 

0.2076 

0.25435 

0.23623 

0.17863 

6 

0.20703 

0.25353 

0.2341 

0.17802 

7 

0.20647 

0.25272 

0.232 

0.17741 

8 

0.20596 

0.25191 

0.22994 

0.17681 

9 

0.20546 

0.25112 

0.22793 

0.17621 

10 

0.20492 

0.25036 

0.22596 

0.17562 

Som. 

0.25189 

0.34105 

0.27768 

0.24917 

Table  XXIV.  Error  norms  (IV.  14)  for  Givoli-Neta  NRBCs,  J  G  1 ...  10,  gravity- 
influenced  system  (IV. 31),  open-air  domain,  left-to-right  advection.  The  error  norm 
from  using  a  Sommerfeld  boundary  condition  is  included  for  comparison. 


but  the  rest  of  the  implementation  is  fully  portable  to  any  other  linear  hrst-order 
system. 

In  the  next  chapter  we  consider  a  third  class  of  NRBC,  one  which  was  recently 
developed  in  the  literature,  which  we  shall  extend  to  our  equation  system. 


112 


VI.  HAGSTROM-WARBURTON  NRBCS 
FOR  THE  LINEARIZED  EULER  EQUATIONS 


A.  INITIAL  IMPLEMENTATION  FOR  FIRST-ORDER 
SYSTEMS 
1.  Derivation 

We  now  consider  a  third  class  of  NRBC,  the  Hagstrom-Warbnrton  anxiliary 
variable  techniqne,  and  apply  it  to  the  linearized  Enler  eqnations.  We  proceed  simi¬ 
larly  to  Sec.  6  of  [57],  which  docnmented  the  concept’s  extension  to  symmetric  hrst- 
order  systems.  While  onr  system  can  be  made  symmetric  [117],  it  is  actnally  not 
necessary  to  do  so,  and  some  nnmerical  experiments  have  shown  that  the  resnlts  are 
either  the  same  as  the  following  straightforward  approach,  worse,  or  nnstable,  depend¬ 
ing  on  the  domain  conhgnration.  The  vector  form  of  (III. 33)  is  a  fairly  straightforward 
extension. 


dt(pi  =  aodt<f  +  dfi0 

ajdt(pj+i  -  dri(pj+i  =  ajdt(pj+i  +  8^0 j  (VI.  1) 

0J+1  =  0 

We  remove  the  wave  speed  c  from  the  eqnations  becanse  the  wave  speeds  will  instead 
be  replaced  by  the  eigenvalnes  of  the  coefficient  matrix  corresponding  to  the  normal 
derivative.  As  with  the  preceding  chapter  for  the  Givoli-Neta  formnlation,  we  begin 
with  a  nsefnl  lemma: 

Lemma  VI. 1  Let  L{(p)  =  dt(p  +  Adx(p  +  Bdytf  —  C(p.  The  auxiliary  variables  ipj 
satisfy 

m)  =  0  (VI.2) 

for  all  j  >  1.  The  only  exception  is  when  both  of  the  following  conditions  apply: 

•  n  =  Lk 

•  ^  0 
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Proof.  Apply  the  differential  operator  L  to  (VI. la),  giving 


L  [dtipi)  =  L  {aodt0  +  dfi0) 

=  (aodt(L{D))+dn{L{D))) 

dt{L{ip,))  =  0  (VI.3) 

Since  0i  is  initially  zero,  L{0i)  is  also  zero,  and  thns  it  is  always  zero.  Proceeding  by 
indnction  we  now  apply  L  to  (VI.  lb)  and  get 

L  {(ijdiip dfi^p L  {o,jdi(pj 

=  ajdt{L{(pj))  +  dri{L{(pj)) 

=  0  (VI.4) 

We  then  integrate  along  the  incoming  characteristic  to  get  L{(pjj^i),  which  mnst  be 
zero  becanse  the  incoming  characteristics  are  zero  □ 

Once  again,  as  with  Lemma  V.l,  the  proof  relies  on  the  fact  that  L{dfi(p)  = 
dfi  [L  {(p)).  On  the  top  and  bottom  of  a  gravity-inffnenced  domain,  that  is  not  trne, 
becanse  the  coefficient  matrices  also  depend  on  We  will  address  this  issne  in  Sec.  C. 

As  with  the  Givoli-Neta  formnlation,  we  wonld  like  to  get  rid  of  the  normal 
derivative  terms  in  (VI.  1).  For  definiteness,  the  following  derivation  is  applied  to  a 
right-hand-side  bonndary,  that  is,  n  =  i.  Simply  make  the  appropriate  changes  for 
the  other  bonndaries  (except  a  vertical  bonndary  in  the  presence  of  gravity,  as  we 
have  noted).  Left-mnltiply  (VI. la)  by  A  and  snbtract  0  =  L{(p)  —  D  to  get 

Adt(pi  =  {ttoA  —  I)  dt(p  —  Bdyif  -f-  Cip  +  D  (VI. 5) 

Similarly,  left-mnitiply  (VI. lb)  by  A  and  add  =  —L{(pj)  =  0  to  get 

(ojA  -f- 1)  dtPj+i  +  BdyPj+i  -  Cpj+i  =  {ajA  -  I)  dtPj  -  BdyPj  +  Cpj  (VI. 6) 

Consider  the  eigenvalne  and  eigenvector  matrices  for  the  Givoli-Neta  NRBCs 
(V.7).  Left-mnitiply  (VI. 5)  and  (VI. 6)  by  R~^  and  make  the  snbstitntion  =  R~^(pj. 
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We  get 


—  (aoA  —  I)  dt^o  —  Bdy^o  +  C^o  +  D 

(a, A  +  I)  dtCj+i  +  -  Cg+i  =  (a, A  -  I)  dtC,  -  Bdyi,  +  Cg  ,  (VI.7) 

where  B,  C  and  D  are  all  defined  by  (V.9).  We  also  insert  R~^R  as  needed  to 
change  the  (pj  terms  to  ^j’s.  We  have  a  self-contained  system  of  eqnations  dehned 
solely  on  the  bonndary.  Hence,  with  the  bonndary  pertnrbations  initially  zero,  these 
eqnations  will  always  remain  zero.  So  we  need  to  incorporate  interior  valnes  into 
onr  formnlation  in  order  to  have  something  more  than  a  rather  convolnted  Dirichlet 
condition.  Following  the  approach  in  [57],  we  replace  the  fonr-eqnation  trnncation 
condition  (pj+i  =  0  with  a  more  characteristic-based  approach.  For  the  ontgoing 
characteristics,  we  nse  the  interior  scheme 

'dtio  +  Aa,go  +  Bdyio  =  Cio  +  d]^  ,  (VI.8) 

where  the  -f-  snbscript  indicates  that  we  only  consider  the  resnlting  eqnations  corre¬ 
sponding  to  the  ontgoing  characteristics.  This  gives  ns  one  or  three  eqnations  (de¬ 
pending  on  the  sign  of  Uq)  for  onr  system  of  nnknowns.  For  the  remaining  nnknowns, 
we  apply  the  trnncation  condition  to  the  incoming  characteristics  only,  i.e., 

[6h-i  =  0]_  (VI.9) 

These  two  eqnations,  in  conjnnction  with  (VI.7),  dehne  a  system  of  4(  J-l-2)  eqnations 
for  4(  J  -|-  2)  nnknowns  ^j.  Collecting  the  time  derivative  terms  onto  one  side  gives  ns 
the  explicit  system: 


dt^o 

A^ig  +  (/  —  aoA)  dtio 
(/  -|-  cijA)  dt^j+i  +  (/  —  cijA)  dt^j 

[gj+i 


-Ad Jo  -  Bdyjo  +  {b  {d,R-^)  R  +  C)^o  +  d]^ 
-Bdyjo  +  {b  (d,R-^)  R  +  C)io  +  D 
-Bdy\J+i  +  (b  {d,Rj  R  +  C)  g+i 
^  -Bdy\J  +  (B{d,RjR  +  C)i, 

0]_  ,  (VI.IO) 
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where  the  subscript  y\z  denotes  a  partial  derivative  with  respect  to  either  y  or  de¬ 
pending  on  the  domain.  The  dzR~^  term  is  non-zero  only  when  gravity  is  considered. 

Note  that  the  presence  of  A  on  the  left-hand  side  of  (VI.  10b)  requires  that  we 
again  have  A  (and  thus  A)  non-singular,  as  with  the  Givoli-Neta  NRBC.  We  solve  the 
problem  here  the  same  way.  Impose  a  false  advection  uo  =  e  ^  0  small  enough  that 
the  system  will  not  drift  across  a  single  grid  point  over  the  duration  of  the  simulation; 
thus, 

St 

ll«oll  <  5—  (VI.ll) 

^^max 

where  Sx  is  the  grid  spacing,  tmax  is  the  duration  of  the  simulation,  and  we  include  a 
factor  of  two  to  ensure  that  there  will  be  no  drift  even  with  round-off  errors. 

2.  Numerical  Examples 

Let  us  consider  a  few  numerical  examples.  We  use  our  usual  inhnite-channel 
example,  considering  the  same  six  cases  as  Sec.  V.3,  always  with  J  G  1 ...  10.  Ta¬ 
ble  XXV  shows  the  state  variable  error  norms  for  the  basic  system  with  no  advection. 
Table  XXVI  shows  the  error  norms  for  the  basic  system  with  left-to-right  advec¬ 
tion.  Table  XXVII  shows  the  error  norms  for  the  Coriolis-influenced  system  with 
no  advection.  Table  XXVIII  shows  the  error  norms  for  the  Coriolis-influenced  sys¬ 
tem  with  left-to-right  {uq  =  100^)  advection.  Table  XXIX  shows  the  error  norms 
for  the  gravity-influenced  system  with  no  advection  (recall  that  we  can  handle  grav¬ 
ity  when  our  open  boundaries  are  on  the  left  and/or  right,  but  not  on  the  top  or 
bottom).  Table  XXX  shows  the  error  norms  for  the  gravity-influenced  system  with 
left-to-right  advection.  For  comparison,  each  table  includes  the  error  norm  resulting 
from  using  a  Sommerfeld  condition.  For  the  “no  advection”  cases  we  use  the  false 
advection  uq  =  1.3889^.  In  every  case,  we  make  the  simplihcation  aj  =  1  for  all  j. 
Hagstrom  et  ah  [53],  citing  an  analysis  by  Diaz  and  Joly  [23],  show  that  this  choice  is 
always  adequate.  We  will  also  make  this  simplihcation  in  all  subsequent  derivations 
for  these  NRBCs.  Fig.  40  shows  the  state  variable  p  at  the  end  of  the  run  for  J  =  10 
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for  the  gravity-influenced  system  with  left-to-right  advection.  Note  again  the  strong 
reflections  which  accompany  the  characteristic-based  implementations  [72], 


Pressure  pefturt>atior^.  reference  solution 
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Figure  40.  The  solution  for  the  pressure  p  using  J  =  10,  Hagstrom-Warburton 
NRBCs,  gravity-influenced  system  (IV. 31),  inhnite  horizontal  channel,  left-to-right 
advection.  (Top)  Reference  solution;  the  area  corresponding  to  the  computed  so¬ 
lution  is  contained  between  the  vertical  lines.  (BL)  Computed  solution  using  NR¬ 
BCs.  (BC)  Reference  solution  domain  corresponding  to  NRBC  solution  domain. 
(BR)  Delta  between  computed  and  reference  solutions,  with  error  norm  computed  by 
(IV.14). 
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J 

E, 

Eu 

Ep 

1 

0.12441 

0.25176 

0.10978 

0.12441 

2 

0.12353 

0.24955 

0.10744 

0.12353 

3 

0.12264 

0.24735 

0.10619 

0.12264 

4 

0.12175 

0.24519 

0.10511 

0.12175 

5 

0.12086 

0.24305 

0.10407 

0.12086 

6 

0.11999 

0.24095 

0.10304 

0.11999 

7 

0.11912 

0.23887 

0.10203 

0.11912 

8 

0.11827 

0.23683 

0.10103 

0.11827 

9 

0.11742 

0.23482 

0.10004 

0.11742 

10 

0.11658 

0.23283 

0.099065 

0.11658 

Som. 

0.21847 

0.53616 

0.12503 

0.21847 

Table  XXV.  Error  norms  (IV.  14)  for  Hagstrom-Warburton  NRBCs  for  basic  system 
(IV.  1)  with  J  G  1 ...  10  in  an  infinite  channel  with  no  advection.  Error  norms  from 
nsing  a  Sommerfeld  radiation  condition  are  inclnded  for  comparison 


J 

E, 

Eu 

Eu 

Ep 

1 

0.13788 

0.2047 

0.14207 

0.13788 

2 

0.1367 

0.20297 

0.13915 

0.1367 

3 

0.13553 

0.20126 

0.13632 

0.13553 

4 

0.13437 

0.19958 

0.13358 

0.13437 

5 

0.13322 

0.19791 

0.13092 

0.13322 

6 

0.13209 

0.19626 

0.12835 

0.13209 

7 

0.13096 

0.19462 

0.12586 

0.13096 

8 

0.12985 

0.19301 

0.12344 

0.12985 

9 

0.12874 

0.19141 

0.1211 

0.12874 

10 

0.12765 

0.18983 

0.11884 

0.12765 

Som. 

0.21817 

0.29414 

0.17727 

0.21817 

Table  XXVI.  Error  norms  (IV.  14)  for  Hagstrom-Warbnrton  NRBCs  for  basic  system 
(IV.  1)  with  J  G  1 ...  10  in  an  infinite  channel  with  left-to-right  advection.  Error 
norms  from  nsing  a  Sommerfeld  radiation  condition  are  inclnded  for  comparison 
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J 

E, 

Eu 

Ey 

Ep 

1 

0.12441 

0.25176 

0.10979 

0.12441 

2 

0.12353 

0.24955 

0.10745 

0.12353 

3 

0.12264 

0.24736 

0.10621 

0.12264 

4 

0.12175 

0.24519 

0.10513 

0.12175 

5 

0.12086 

0.24306 

0.10408 

0.12086 

6 

0.11999 

0.24095 

0.10306 

0.11999 

7 

0.11912 

0.23888 

0.10204 

0.11912 

8 

0.11827 

0.23684 

0.10104 

0.11827 

9 

0.11742 

0.23482 

0.10005 

0.11742 

10 

0.11658 

0.23284 

0.09908 

0.11658 

Som. 

0.21847 

0.53616 

0.12503 

0.21847 

Table  XXVII.  Error  norms  (IV.  14)  for  Hagstrom-Warburton  NRBCs  for  Coriolis- 
influenced  system  (IV.  15)  with  J  G  1 ...  10  in  an  infinite  channel  with  no  advection. 
Error  norms  from  using  a  Sommerfeld  radiation  condition  are  included  for  comparison 


J 

E, 

Eu 

Ey 

Ep 

1 

0.13855 

0.31144 

0.21862 

0.13855 

2 

0.13753 

0.30948 

0.21789 

0.13753 

3 

0.13654 

0.3076 

0.21725 

0.13653 

4 

0.13557 

0.3058 

0.21668 

0.13557 

5 

0.13464 

0.30407 

0.21617 

0.13464 

6 

0.13372 

0.3024 

0.21573 

0.13372 

7 

0.13284 

0.3008 

0.21534 

0.13284 

8 

0.13198 

0.29926 

0.215 

0.13197 

9 

0.13113 

0.29777 

0.21471 

0.13113 

10 

0.13032 

0.29633 

0.21446 

0.13032 

Som. 

0.1929 

0.32104 

0.15953 

0.1929 

Table  XXVIII.  Error  norms  (IV.  14)  for  Hagstrom-Warburton  NRBCs  for  Coriolis- 
influenced  system  (IV.  15)  with  J  G  1 ...  10  in  an  inhnite  channel  with  left-to-right 
advection.  Error  norms  from  using  a  Sommerfeld  radiation  condition  are  included  for 
comparison 
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J 

E, 

Eu 

Ev 

Ep 

1 

0.12896 

0.25097 

0.11173 

0.12835 

2 

0.12807 

0.24879 

0.10934 

0.12745 

3 

0.12714 

0.24663 

0.10807 

0.12652 

4 

0.12622 

0.24449 

0.10697 

0.1256 

5 

0.1253 

0.24238 

0.10591 

0.12469 

6 

0.12439 

0.24031 

0.10487 

0.12379 

7 

0.12349 

0.23826 

0.10384 

0.12289 

8 

0.12259 

0.23624 

0.10283 

0.12201 

9 

0.12171 

0.23425 

0.10183 

0.12113 

10 

0.12084 

0.2323 

0.10084 

0.12026 

Som. 

0.22442 

0.52603 

0.12775 

0.22374 

Table  XXIX.  Error  norms  (IV.  14)  for  Hagstrom-Warburton  NRBCs  for  gravity- 
influenced  system  (IV. 31)  with  J  G  1 ...  10  in  an  infinite  channel  with  no  advection. 
Error  norms  from  using  a  Sommerfeld  radiation  condition  are  included  for  comparison 


J 

E, 

Eu 

Ey 

Ep 

1 

0.13899 

0.20395 

0.14268 

0.13852 

2 

0.13779 

0.20223 

0.13977 

0.13733 

3 

0.1366 

0.20053 

0.13696 

0.13615 

4 

0.13543 

0.19885 

0.13423 

0.13498 

5 

0.13426 

0.19718 

0.13158 

0.13382 

6 

0.13311 

0.19554 

0.12902 

0.13267 

7 

0.13197 

0.19391 

0.12654 

0.13154 

8 

0.13084 

0.1923 

0.12414 

0.13041 

9 

0.12972 

0.19071 

0.12181 

0.1293 

10 

0.12861 

0.18913 

0.11956 

0.1282 

Som. 

0.21869 

0.28969 

0.17855 

0.21802 

Table  XXX.  Error  norms  (IV.  14)  for  for  Hagstrom-Warburton  NRBCs  for  gravity- 
influenced  system  (IV.31)  with  J  G  1 ...  10  in  an  inhnite  channel  with  left-to-right 
advection.  Error  norms  from  using  a  Sommerfeld  radiation  condition  are  included  for 
comparison 
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B.  CORNER  CONDITIONS 
1.  Derivation 

We  can  use  the  same  methods  as  in  Sec.  V.l  to  apply  these  NRBCs  to  the 
corners  of  an  open  domain,  in  the  absence  of  gravity.  Using  the  same  notation  as 
before,  our  characteristic-based  NRBC  scheme  on  three  open  sides  of  a  half-plane  is 
given  by 


dtC+  =  -AdnC+  -  Tdr0+ 
Adti,  +  {I-A)dti  =  -fdr0 
(/  +  A)aig+i  +  (/-A)aig  = 

ij+i-  =  0 


where  we  use 


dn 

Left  :  <  dr 

T 

.. 

dn 

Right  :  <  dr 
T 

.. 

dn 

Top  :  <  dr 

T 

dn 

Top-left  :  <  dr 

T 

dn 

Top-right  :  dr 

T 


-dx 

dz 

B 

dx 

dz 

B 

dz 

dx 

A 


(—dx  +  dz)  / -\/2 
(dx  +  dz)  / 

(A  +  B)/V2 

(dx  +  dz)  / 

(dx  -  dz)  /v^ 
(A  -  B)  7^2 


(VL12) 


(VL13) 
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2.  Numerical  Examples 

Let  us  now  consider  some  examples.  Since  we  cannot  include  the  effects  of 
gravity  on  the  top  of  the  domain,  we  consider  only  four  cases:  The  basic  system 
and  the  system  influenced  by  Coriolis  forces,  with  and  without  horizontal  advec- 
tion.  We  use  the  same  false  advection  as  in  Sec.  2  (horizontal)  and  V.2  (vertical). 
Table  XXXI  shows  the  error  norms  for  J  G  1 ...  10  for  the  basic  system  with  no  ad¬ 
vection.  Table  XXXII  shows  the  error  norms  for  J  G  1 ...  10  for  the  basic  system  with 
left-to-right  advection.  Table  XXXIII  shows  the  error  norms  for  J  G  1 ...  10  for  the 
Coriolis-influenced  system  with  no  advection.  Table  XXXIV  shows  the  error  norms 
for  J  G  1 ...  10  for  the  Coriolis-influenced  system  with  left-to-right  advection.  Fig.  41 
shows  an  example  of  the  horizontal  velocity  perturbation  u  for  the  basic  no-advection 
test.  Note  the  reflections  are  stronger  at  the  bottom  corners  than  the  top.  These 
reflections  are  the  result  of  waves  reflecting  off  the  hard  wall  on  the  bottom  before 
impacting  the  left  and  right  boundaries;  thus,  their  incidence  angles  are  greater  than 
45°.  At  the  top  corners,  these  waves  pass  through  the  open  top  with  the  smaller  inci¬ 
dence  angle,  so  there  is  less  reflection  evident.  If  we  plot  the  values  of  Table  XXXIV, 
we  see  that  the  error  norms  decrease  exponentially  as  J  increases  by  2  (Fig.  42,  but 
note  that  v  increases  with  J  >  7).  The  reason  for  the  saw-tooth  pattern  is  unknown, 
but  it  is  only  present  in  the  Coriolis-l-advection  case;  the  other  cases  show  exponential 
reductions  with  each  increase  in  J. 
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X-direclion  vel^ity  perturtation.  reference  solution 


-50  0  so  100  ISO  200 


X-direction  velocity  perturbation  Truncated  reference  solution  Solution  delta 
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Error  nofm-0-2758 
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_ 
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0.02 
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-0.04 


Figure  41.  The  solution  for  the  horizontal  velocity  u  using  J  =  10,  Hagstrom- 
Warburton  NRBCs,  basic  system  (IV. 1),  open  half-plane,  no  advection.  (Top)  Ref¬ 
erence  solution;  the  area  corresponding  to  the  computed  solution  is  contained  in  the 
bottom-center  box.  (BL)  Computed  solution  using  NRBCs.  (BC)  Reference  solution 
domain  corresponding  to  NRBC  solution  domain.  (BR)  Delta  between  computed  and 
reference  solutions,  with  error  norm  computed  by  (IV.  14). 


J 

E, 

Eu 

Ey 

Ep 

1 

0.17603 

0.30011 

0.13089 

0.17604 

2 

0.17485 

0.29724 

0.12848 

0.17485 

3 

0.17364 

0.29438 

0.12699 

0.17365 

4 

0.17244 

0.29157 

0.12565 

0.17245 

5 

0.17125 

0.28881 

0.12435 

0.17126 

6 

0.17007 

0.2861 

0.12307 

0.17007 

7 

0.1689 

0.28345 

0.12182 

0.1689 

8 

0.16773 

0.28085 

0.12059 

0.16774 

9 

0.16658 

0.2783 

0.11938 

0.16658 

10 

0.16543 

0.2758 

0.11819 

0.16543 

Som. 

0.33324 

0.71942 

0.2414 

0.33325 

Table  XXXI.  Error  norms  (IV.  14)  for  Hagstrom-Warburton  NRBCs  for  basic  system 
(IV.  1)  with  J  G  1 ...  10  in  an  open  half-plane  with  no  advection.  Error  norms  from 
using  a  Sommerfeld  radiation  condition  are  included  for  comparison 
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J 

E, 

Eu 

Ey 

Ep 

1 

0.16823 

0.20592 

0.18956 

0.16823 

2 

0.16661 

0.20403 

0.18543 

0.16661 

3 

0.16501 

0.20217 

0.18143 

0.16501 

4 

0.16342 

0.20034 

0.17755 

0.16342 

5 

0.16186 

0.19854 

0.1738 

0.16186 

6 

0.16031 

0.19675 

0.17016 

0.16031 

7 

0.15878 

0.19499 

0.16664 

0.15878 

8 

0.15727 

0.19326 

0.16323 

0.15727 

9 

0.15577 

0.19154 

0.15993 

0.15577 

10 

0.15429 

0.18985 

0.15673 

0.1543 

Som. 

0.2675 

0.32527 

0.2411 

0.2675 

Table  XXXII.  Error  norms  (IV.  14)  for  Hagstrom-Warburton  NRBCs  for  basic  system 
(IV.  1)  with  J  G  1 ...  10  in  an  open  half-plane  with  left-to-right  advection.  Error 
norms  from  nsing  a  Sommerfeld  radiation  condition  are  inclnded  for  comparison 


J 

E, 

Eu 

Eu 

Ep 

1 

0.17604 

0.30024 

0.13095 

0.17604 

2 

0.17484 

0.29737 

0.12852 

0.17484 

3 

0.17364 

0.29451 

0.12704 

0.17365 

4 

0.17243 

0.29169 

0.12569 

0.17243 

5 

0.17125 

0.28893 

0.12439 

0.17126 

6 

0.17006 

0.28622 

0.12311 

0.17006 

7 

0.16889 

0.28356 

0.12186 

0.1689 

8 

0.16772 

0.28096 

0.12063 

0.16772 

9 

0.16657 

0.2784 

0.11943 

0.16658 

10 

0.16541 

0.27591 

0.11824 

0.16542 

Som. 

0.33324 

0.71943 

0.2414 

0.33325 

Table  XXXIII.  Error  norms  (IV.  14)  for  Hagstrom-Warbnrton  NRBCs  for  Coriolis- 
inflnenced  system  (IV.  15)  with  J  G  1 ...  10  in  an  open  half-plane  with  no  advection. 
Error  norms  from  nsing  a  Sommerfeld  radiation  condition  are  inclnded  for  comparison 
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J 

E, 

Eu 

Ey 

Ep 

1 

0.20445 

0.35374 

0.261 

0.20444 

2 

0.20472 

0.35221 

0.26198 

0.20471 

3 

0.20277 

0.34939 

0.26073 

0.20276 

4 

0.20307 

0.34806 

0.26177 

0.20306 

5 

0.20117 

0.34542 

0.26059 

0.20116 

6 

0.20152 

0.34424 

0.26168 

0.2015 

7 

0.19966 

0.34175 

0.26057 

0.19965 

8 

0.20003 

0.34072 

0.26169 

0.20002 

9 

0.19822 

0.33836 

0.26064 

0.19821 

10 

0.19862 

0.33744 

0.26179 

0.19861 

Som. 

0.26015 

0.34856 

0.20456 

0.26014 

Table  XXXIV.  Error  norms  (IV.  14)  for  Hagstrom-Warburton  NRBCs  for  Coriolis- 
influenced  system  (IV.  15)  with  J  G  1 ...  10  in  an  open  half-plane  with  left-to-right 
advection.  Error  norms  from  using  a  Sommerfeld  radiation  condition  are  included  for 
comparison 


Density  perturbation  error  norm  X-velocity  perturbation  error  norm 


Y-velocity  perturbation  error  norm  Pressure  perturbation  error  norm 


Figure  42.  Logarithmic  plot  of  error  norms  (IV.  14)  for  Hagstrom-Warburton  NRBC, 
J  G  2  ...  10,  open  half-plane,  Coriolis,  left-to-right  advection.  (TL)  Error  norm  for  p. 
(TR)  Error  norm  for  u.  (BL)  Error  norm  for  v.  (BR)  Error  norm  for  p. 
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C.  INCOMPATIBILITY  WITH  GRAVITY 


Problems  arise,  however,  when  we  try  to  incorporate  the  effects  of  gravity  on 
the  top  bonndary.  (We  can  inclnde  it  on  P/,  and  Pr  withont  any  difficnlty.)  Using 
the  same  approach  as  for  the  Givoli-Neta  NRBCs,  we  apply  the  differential  operator 
L  to  both  sides  of  (VI. la),  beginning  the  process  of  removing  the  normal  derivatives: 


L  (5t(pi)  =  L  (ao5t(p  + 


aodt{L{(p))  +  L{d^(p) 

aodtD  +  dt^(f  +  -  Cd^(p 

/ 

dz  {dt(p  +  Ad^ip  +  Bdz0  -  C(p) 

-  {dzA)  (p  -  {dzB)  (p  +  {dzC)  (p 
d,  (L  (<p))  -  {3^)0-  {d,B)0+  {d,C)  (p  .  (VI.  14) 


Even  with  the  simplifying  assnmption  of  an  exponentially-decaying  atmosphere,  we 
still  have 


dt{,L  [tpi)) 


dt{L  {0i)) 


dzD  +  aAzff  +  aBz0  —  aC(p 

-aD  +  a  (L  {tp)  -  dt0  -  Uod^0  -  Wodz<p) 

-aD  +  aD  -  a  {dt(p  +  uodx(p  +  wodz(p) 

-adF(p  .  (VI.  15) 


To  continne  this  derivation,  we  mnst  integrate  the  Lagrangian  flow  derivative  of  (p  over 
time.  The  resnlt  is  then  applied  to  find  L  (<^2),  with  increasingly  convolnted  combina¬ 
tions  of  Lagrangian  derivatives  and  time  integrals.  This  approach  is  not  satisfactory. 
Hence,  we  reqnire  a  different  approach  if  we  wish  to  incorporate  gravitational  effects 
into  the  Hagstrom-Warbnrton  NRBCs.  We  have  not  fonnd  snch  an  approach,  bnt 
fntnre  research  may  reveal  a  snitable  techniqne. 
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D.  WAVE-LIKE  IMPLEMENTATION 
1.  Derivation 

The  results  from  this  implementation  are  not  bad,  but  the  improvement  with 
larger  J  is  slight.  Furthermore,  the  reflections  downstream  of  the  advection  are  trou¬ 
blesome.  Perhaps  we  can  do  better.  Looking  at  (VI.  1),  we  see  a  hint  of  the  problem. 
Although  Hagstrom  and  Warburton  do  not  provide  a  physical  interpretation  of  their 
auxiliary  variable  formulation,  inspection  of  (VI. 1)  shows  that  we  can  conceive  of  the 
second  equation  in  the  following  manner:  The  outgoing  characteristic  of  (pj  is  in  some 
sense  paired  with  the  incoming  characteristic  of  Pj+i-  However,  when  we  use  these 
characteristic  variables  for  a  hrst-order  system,  we  contradict  this  interpretation  of 
the  auxiliary  variables,  since  each  variable  has  only  one  characteristic  (either  incom¬ 
ing  or  outgoing),  not  two.  If  we  could  contrive  a  second  characteristic  for  each  state 
variable,  it  might  improve  our  results.  Let  us  instead  apply  the  NRBC  to  each  state 
variable  individually.  Instead  of  Lemma  VI.l,  we  use  the  following: 


Lemma  VI. 2  Each  state  variable  p  G  {p,  has  a  solution  which  satisfies  the 

acoustic  wave  equation 


dttp  =  VV  , 


(VL16) 


where 


Co  +  Uo 

on  Tg 

Co  ~ 

onTw 

Co  +  Vo 

on  Tat 

o 

1 

O 

on  T^ 

(VI.  17) 


The  derivation  is  given  in  Appendix  C.  With  this  fact,  we  can  then  implement  the 


auxiliary  variables  exactly  as  given  in  [57],  taking  each  state  variable  individually. 
Upon  removing  the  normal  derivatives,  we  are  left  with  the  following  system  of  equa¬ 


tions: 


—a^dtp  +  dtPi 
2ai(l  -  al)dttp 
-|-(1  +  a‘1  +■  2aoai)dttPi  ^ 
+  (1  ~  (^\)dttP2 


=  CpdxP 

/ 

cl{2aidyyP 
+  dyyPl  +  dyyP2) 
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aj{l- 

-\-{aj-i  +  +  (ij-i(ij)dtt‘fj 

+%-i(l  —  aj)du‘fj+i 

ap{l-  al,_^)dtt(pp-i 
+  {ap-i  +  cip){l  +  ap-iap)dtt<fp 


>  = 


>  = 


cl{aAy^j-i 

+  (aj-i  +  aj)dyy(pj 
-\-aj-idyy(pjpi) 

cliapdyyifp-i 
+  (ap_i  +  ap)dyy(pp) 


(VL18) 


where  =  Co  +  Uq.  When  we  make  the  simplihcation  a,  =  1,  this  system  reduces  to 


-dtif  +  dtifi 

4:du‘fi 

4^dtt(Pj 

^dtt<fp 


Cpdx(p 

/ 

Cl{‘^dyyip 

< 

+  dyyipi  +  dyy(f2) 

/ 

+2dyy(pj 

+  dyy(pjpl) 

V 

/ 

cl{dyy^P-l 

< 

-\-2dyy(pp')  , 


(VL19) 


This  system  is  the  NRBC  on  Tp.  On  Tw  replace  dx^  in  (VI. 19a)  with  —dx^p,  and 
use  Cxu  =  Co  —  Uq- 


2.  Numerical  Examples 

If  we  run  the  same  infinite-channel  example  as  in  Sec.  1.2,  using  the  wave-based 
auxiliary  variable  approach  from  the  preceding  section,  we  get  the  results  illustrated  in 
Fig.  43  with  the  error  norms  for  all  state  variables  given  in  Table  XXXV.  If  we  use  the 
non-zero  mean  flow  uq  =  100^,  we  get  the  error  norms  in  Table  XXXVI.  Comparing 
Tables  XXXV  and  XXXVI  to  Tables  XXV  and  XXVI,  respectively,  we  see  that  this 
new  version’s  error  norms  are  approximately  half  the  old  version’s.  We  note  also  that 
the  new  method’s  error  norms  show  almost  no  improvement  for  J  >  5.  This  “error 
floor”  has  been  observed  in  other  auxiliary  variable  implementations  for  scalar  wave 
equations  [36].  Although  the  old  method  shows  continual  (albeit  slow)  improvement, 
the  error  norms  do  not  match  those  of  the  new  method  until  J  ~  500,  which  runs  far 
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Pressure  perturtation.  reference  solution 
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Figure  43.  The  solution  for  the  pressure  p  using  J  =  10,  wave-like  Hagstrom- 
Warburton  NRBCs  in  an  inhnite  channel,  basic  system  with  no  advection.  (Top)  Ref¬ 
erence  solution;  the  area  corresponding  to  the  computed  solution  is  contained  between 
the  vertical  lines.  (BL)  Computed  solution  using  NRBCs.  (BC)  Reference  solution 
domain  corresponding  to  NRBC  solution  domain.  (BR)  Delta  between  computed  and 
reference  solutions,  with  error  norm  computed  by  (IV.  14). 

slower  than  the  new  method’s  J  =  5.  While  this  result  is  promising,  experiments  with 
longer  integration  times  show  that  this  wave-like  implementation  is  less  stable  than 
the  characteristic-based  method;  hence,  we  do  not  pursue  its  development  further. 


E.  SUMMARY 

This  chapter  concludes  our  development  of  NRBCs  for  the  linearized  Euler 
equations.  For  each  NRBC,  we  have  derived  its  implementation  and  demonstrated 
its  effectiveness  for  short  time-integrations.  In  the  next  two  chapters,  we  will  consider 
the  question  of  longer  time-integrations  and  the  relative  strengths  and  weaknesses  of 
each  method. 
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J 

E, 

Eu 

Ep 

1 

0.059866 

0.16273 

0.038642 

0.059866 

2 

0.059442 

0.16184 

0.038079 

0.059443 

3 

0.059306 

0.16156 

0.037884 

0.059306 

4 

0.059252 

0.16144 

0.037804 

0.059252 

5 

0.059236 

0.16139 

0.037775 

0.059237 

6 

0.059234 

0.16138 

0.037765 

0.059235 

7 

0.059235 

0.16138 

0.037761 

0.059236 

8 

0.059236 

0.16138 

0.03776 

0.059237 

9 

0.059237 

0.16138 

0.037759 

0.059237 

10 

0.059237 

0.16138 

0.037759 

0.059238 

Table  XXXV.  Error  norms  (IV.  14)  for  J  G  1 ...  10,  wave-like  Hagstrom-Warburton 
NRBCs,  infinite  channel,  basic  system  with  no  advection 


J 

E, 

Eu 

Eu 

Ep 

1 

0.053965 

0.078277 

0.040606 

0.053965 

2 

0.053561 

0.077765 

0.040013 

0.053561 

3 

0.05347 

0.077625 

0.039873 

0.05347 

4 

0.053457 

0.077589 

0.039842 

0.053457 

5 

0.053456 

0.077575 

0.039836 

0.053456 

6 

0.053455 

0.077568 

0.039835 

0.053455 

7 

0.053454 

0.077564 

0.039835 

0.053454 

8 

0.053453 

0.077562 

0.039835 

0.053453 

9 

0.053453 

0.077561 

0.039834 

0.053453 

10 

0.053453 

0.077561 

0.039834 

0.053453 

Table  XXXVI.  Error  norms  (IV.  14)  for  J  G  1 ...  10,  wave-like  Hagstrom-Warbnrton 
NRBCs,  inhnite  channel,  basic  system  with  left-to-right  advection 
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VII. 


LONG-TIME  STABILITY 


A.  OBSERVATIONS 

In  this  chapter  we  address  the  stability  of  the  NRBCs  for  long  time-integrations. 
This  issue  is  critical  to  numerical  weather  prediction  and  other  applications  which  re¬ 
quire  time  frames  beyond  that  of  a  single  wave’s  propagation  through  the  domain. 

Higdon  [68]  discusses  the  numerical  stability  of  his  NRBC  formulation,  us¬ 
ing  criteria  developed  by  Kreiss  [81]  and  Gustafsson  et  ah  [52]  (and  Higdon’s  own 
characteristic-based  interpretation  thereof  [63]).  While  his  scheme  meets  the  de- 
hned  stability  criteria,  there  is  still  concern  that  a  scheme  which  is  stable  for  the 
single-variable  Klein-Gordon  equation  is  also  stable  long-term  for  the  equivalent 
linearized  Euler  system.  Long-term  stability,  surprisingly,  has  not  been  much  dis¬ 
cussed  or  demonstrated  in  NRBG  development.  In  the  Givoli-Neta  papers  explor¬ 
ing  the  automated  Higdon  NRBGs  and  the  Givoli-Neta  auxiliary  variable  NRBGs 
[39,  40,  41,  42,  43,  90,  113,  114,  115,  116],  only  [90]  discusses  the  long-time  stability 
of  the  solution,  showing  that  the  automated  Higdon  NRBG  of  order  J  =  10  is  stable 
for  long  time-integrations.  In  the  papers  exploring  the  Hagstrom-Warburton  auxil¬ 
iary  variable  scheme  [57,  37,  53,  55],  long-time  stability  was  claimed  in  [57]  (although 
the  numerical  example  plotted  showed  an  error  norm  that  was  increasing  over  time), 
and  [55]  uses  an  evanescent  mode  correction  (a  second  set  of  auxiliary  variables)  to 
ensure  long-time  stability. 

As  for  the  papers  which  discuss  various  NRBG  techniques  for  the  linearized 
Euler  equations,  only  the  absorbing  layer  methods  [1,  17,  72,  73,  88]  discuss  or  demon¬ 
strate  long-time  stability.  Specihcally,  the  numerical  examples  in  [17]  are  carried  out 
for  long  time-integrations,  and  the  PMLs  [1,  72,  73]  require  additional  terms  to  make 
them  stable  ([88]  claims  medium-term  stability  and  implies  stability  over  long  time- 
integrations). 

All  three  types  of  NRBGs  presented  in  this  dissertation  suffer  from  some  form 
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of  long-time  instability.  One  example  may  be  seen  in  Fig.  44,  which  is  the  state 
variable  u  from  the  nnmerical  example  of  Sec.  IV. B. 3,  with  J  =  10,  rnn  nntil  t  =  100  s. 
As  we  can  see,  the  valnes  on  the  bonndary  have  become  catastrophically  large,  and 
they  have  pollnted  the  domain,  completely  overwhelming  the  actnal  wave  (faintly 
visible  near  the  center  of  the  domain).  Tables  XXXVII-XXXIX  list  the  maximnm 
order  J  for  which  the  given  NRBC  is  stable  for  the  varions  domain  conhgnrations 
for  short  (24  s),  medium  (100  s),  and  long  (10,000  s)  time- integrations.  Table  entries 
marked  “No”  indicate  that  the  conhgnration  was  nnstable  even  for  J  =  1.  The 
npper  limit  to  J  for  the  short-term  stability  mirrors  a  condition  seen  while  prodncing 
the  resnlts  given  in  [90] ;  even  thongh  the  Higdon  scheme  is  theoretically  stable  in  the 
discrete  case,  ronnd-off  errors  in  the  hnite-precision  implementation  lead  to  instability. 
This  same  problem  afflicted  the  examples  given  in  Chapter  IV:  If  the  example  was 
performed  on  a  50  x  50  grid,  it  was  stable  only  np  to  J  =  5;  halving  the  grid  spacings 
enabled  the  stable  resnlts  np  to  J  =  10. 


X'direction  Vetocity  pefturtatron 


Fignre  44.  Plot  of  nnstable  u  in  basic  system  (IV. 1)  with  Higdon  NRBC  J  =  10  in 
a  semi-inhnite  channel  integrated  np  to  t  =  100  s.  Notice  the  faint  wave  crests  near 
the  middle  of  the  domain 
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Domain 

Coriolis 

Short 

Medium 

Long 

shape 

or  gravity 

Advection 

(24  s) 

(100  s) 

(10,000  s) 

Bucket 

None 

None 

10 

8 

4 

Coriolis 

None 

10 

8 

4 

Gravity 

None 

10 

8 

No 

Channel 

None 

L-R 

10 

4 

No 

Coriolis 

L-R 

2 

2 

No 

Gravity 

None 

10 

8 

4 

L-R 

10 

4 

No 

Open 

None 

None 

10 

5 

2 

L-R 

10 

4 

No 

BL-TR 

10 

4 

No 

Coriolis 

None 

10 

5 

2 

L-R 

2 

No 

No 

BL-TR 

2 

No 

No 

Ground 

Gravity 

None 

10 

5 

No 

L-R 

10 

4 

No 

Table  XXXVII.  Higdon  NRBCs,  maximum  stable  order  J  for  various  domain  config¬ 
urations  and  simulation  durations 


For  the  Hagstrom-Warburton  NRBCs,  we  found  short-term  stability  as  high 
as  J  =  40.  While  the  NRBC  still  exhibits  exponential  improvement  (see  Fig.  45), 
the  improvement  is  slight  enough  that  the  computational  overhead  for  such  high  J  is 
unjustihed;  hence,  we  did  not  seek  an  actual  “maximum”  stable  J  but  simply  noted 
that  the  maximum  is  at  least  J  =  40. 
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Domain 

Coriolis 

Short 

Medium 

Long 

shape 

or  gravity 

Advection 

(24  s) 

(100  s) 

(10,000  s) 

Bucket 

None 

None 

15 

7 

No 

Coriolis 

None 

15 

7 

No 

Gravity 

None 

13 

6 

No 

Channel 

None 

L-R 

12 

7 

4 

Coriolis 

L-R 

11 

7 

1 

Gravity 

None 

15 

7 

No 

L-R 

12 

7 

3 

Ground 

None 

None 

15 

7 

No 

L-R 

12 

3 

No 

Coriolis 

None 

15 

7 

No 

L-R 

10 

3 

No 

Gravity 

None 

13 

No 

No 

L-R 

11 

3 

No 

Table  XXXVIII.  Givoli-Neta  NRBCs,  maximum  stable  order  J  for  various  domain 
configurations  and  simulation  durations 


Domain 

Coriolis 

Short 

Medium 

Long 

shape 

or  gravity 

Advection 

(24  s) 

(100  s) 

(10,000  s) 

Bucket 

None 

None 

40+ 

40+ 

No 

Coriolis 

None 

40+ 

40+ 

No 

Channel 

None 

L-R 

40+ 

40+ 

40+ 

Coriolis 

L-R 

40+ 

40+ 

40+ 

Gravity 

None 

40+ 

40+ 

No 

L-R 

40+ 

40+ 

5 

Ground 

None 

None 

40+ 

28 

No 

L-R 

40+ 

3 

No 

Coriolis 

None 

40+ 

27 

No 

L-R 

40+ 

4 

No 

Table  XXXIX.  Hagstrom-Warburton  NRBCs,  maximum  stable  order  J  for  various 
domain  configurations  and  simulation  durations 
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Density  perturbation  error  norm 


X-velocity  perturbation  error  norm 


Y-velocity  perturbation  error  norm  Pressure  perturbation  error  norm 


Figure  45.  Logarithmic  plot  of  error  norms  (IV.  14)  for  J  G  1 . . .  40,  Hagstrom- 
Warburton  NRBCs,  open  half-plane  with  Coriolis,  no  advection.  (TL)  Error  norms 
for  p.  (TR)  Error  norms  for  u.  (BL)  Error  norms  for  v.  (BR)  Error  norms  for  p. 
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B.  SPECULATIONS 

Based  on  observations  from  these  and  other  experiments  dnring  this  research, 
we  can  tentatively  identify  several  possible  instability  sonrces.  These  sonrces,  thongh 
cnrrently  mere  conjectnre  and  specnlation,  will  be  explored  in  fntnre  research. 

1.  Round-off  errors.  Early  tests  of  these  NRBCs  nsed  a  smaller  domain,  only  1 
km  sqnare,  and  the  scheme  was  stable  np  to  J  ~  8.  Increasing  the  domain 
size  and  grid  spacing  by  an  order  of  magnitnde  made  the  scheme  stable  np  to 
J  =  10. 

2.  Normal  derivatives  reaching  too  far  into  the  domain.  The  hrst  nnmerical  ex¬ 
amples,  nsing  the  Higdon  scheme,  were  rnn  on  a  50  x  50  grid.  The  scheme  was 
nnstable  for  J  >  5.  Changing  the  grid  to  100  x  100  made  the  scheme  stable 
np  to  J  =  10.  As  a  resnlt,  the  normal  derivatives  of  the  Higdon  scheme  did 
not  reach  as  far  into  the  domain.  Van  Joolen  [113]  noted  this  same  behavior 
for  the  Klein-Gordon  eqnation. 

3.  Dependency  between  interior  and  NRBC  near  corners.  Althongh  no  points 
depend  on  the  corner  valnes  in  the  Higdon  NRBC  scheme  (see  Fig.  13),  the 
interior  point  closest  to  the  corner  depends  on  two  bonndary  points,  rather 
than  jnst  one.  It  conld  be  that  this  dependence  creates  a  slight  error,  jnst 
enongh  that  it  grows  over  time  and  destabilizes  the  system. 

4.  Frequency  mis-match  between  interior  and  Higdon  NRBC  schemes.  A  reviewer 
for  [21]  noted  that  the  Higdon  NRBC  discretization  scheme  (IV. 10)  cannot  re¬ 
solve  the  shortest  wavelengths  resolvable  by  the  interior  scheme  (IV.  7).  When 
a  wave  strikes  the  open  bonndary,  the  NRBC  compntes  the  valne  based  on 
every  other  point  at  every  other  time  step,  as  if  for  a  wave  with  twice  the 
wavelength.  The  skipped  points  of  the  wave  are  then  nsed  at  the  next  time 
step  to  compnte  that  next  time  step’s  bonndary  valne.  Hence,  a  wave  striking 
the  bonndary  with  wavelength  n  is  resolved  by  the  NRBC  as  a  two-stage  com¬ 
position  of  two  waves,  each  with  wavelength  2z/.  A  stability  analysis  by  one  of 
the  co-anthors  of  [21]  is  in  progress  to  determine  the  impact  of  this  resolntion 
discrepancy. 

5.  Numerically-singular  matrix  computations.  For  the  Givoli-Neta  NRBCs  with 
gravity,  in  either  the  “bncket”  or  open-air  domains,  Matlab  issned  “Matrix  is 
close  to  singnlar  or  badly  scaled”  warnings  for  large  J  and  small  St.  It  appears 
the  left-hand  matrix  nsed  to  solve  the  anxiliary  variables  is  prone  to  nnmeric 
instability. 

6.  Exponentially- growing  solutions.  Althongh  Higdon  [68]  proves  the  stability 
of  his  system  by  demonstrating  that  exponentially-growing  solntions  are  not 
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permitted  for  the  discrete  system  and  boundary  condition,  it  is  possible  that 
such  solutions  are  appearing  in  this  system.  Look  at  Fig.  46,  which  shows  the 
results  of  an  unstable  conhguration  run  until  t  =  100s.  For  each  unstable  J,  the 
oo-norm  of  the  solution  grows  exponentially  after  a  certain  point  in  time.  If  we 
were  to  extend  these  lines  backward,  they  would  intersect  at  approximately 
t  =  0.  Hence,  it  appears  that  there  is  an  initially  small,  but  exponentially 
growing,  solution  within  the  scheme.  Over  time,  this  exponentially  growing 
solution  overwhelms  the  other  interior  values  and  leads  to  the  unstable  system. 

7.  Discrete  reflection  coefficient  greater  than  unity.  Our  analysis  of  the  reflection 
coefficients  (III.  15)  and  (III. 34)  is  based  on  the  partial  derivatives  in  the  con¬ 
tinuous  case.  However,  we  implement  the  system  using  hnite  differences  on 
discrete  variables.  It  is  possible  that  the  “reflection  coefficient”  in  the  discrete 
case  is  dependent  on  more  factors  than  just  the  wave  speed  and  the  NRBC 
estimated  wave  speed,  factors  such  as  the  grid  spacing  and  the  time  step  size. 
Some  conhgurations  may  lead  to  ||i?||  >  1  for  some  waves.  See  Appendix  D 
for  some  initial  analysis.  A  full  analysis  of  this  discrete  reflection  coefficient 
will  be  performed  in  future  research. 


Y-velocity  infinity  norm  over  time 


Pressure  infinity  norm  over  time 


Figure  46.  Logarithmic  plot  of  state  variable  oo-norms  for  J  G  1 ...  11,  Givoli-Neta 
NRBCs,  open  half-plane,  left-to-right  advection,  no  gravity  or  Coriolis,  run  until 
t  =  100s.  (TL)  (X)-norms  for  p.  (TR)  oo-norms  for  u.  (BL)  oo-norms  for  w.  (BR) 
oo-norms  for  p. 


One  result  of  the  Hagstrom-Warburton  table  was  surprising:  The  “bucket” 
configuration  was  always  unstable  long-term,  but  the  infinite  channel  configuration 
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showed  long-term  stability.  This  result  was  not  a  programming  error.  Running  the 
simulation  on  a  horizontal  semi-inhnite  channel  showed  the  same  instability  as  the 
vertical  semi-inhnite  channel.  The  stability  of  the  boundary  condition  was  impacted 
by  the  type  of  boundary  on  the  opposite  side.  Until  a  formal  analysis  is  undertaken 
to  explain  this  result,  we  offer  here  a  tentative  hypothesis.  The  instability  comes 
from  waves  whose  rehection  coefficients  are  slightly  greater  than  unity  when  striking 
an  open  boundary  on  one  side.  On  an  open  boundary  on  the  other  side,  the  wave 
speed  and  NRBC  wave  speed  combine  to  cause  a  rehection  coefficient  less  than  unity; 
however,  if  the  other  side  is  a  hard  wall,  then  its  rehection  coefficient  is  exactly  unity, 
and  so  the  undiminished  wave  is  totally  rehected  back  to  the  open  boundary,  where 
it  is  again  rehected  and  magnihed  back  toward  the  hard  wall,  and  so  forth,  growing 
in  magnitude  with  each  pass. 
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VIII. 


SUMMARY  AND  COMPARISONS 


Having  developed  three  distinct  NRBC  methods,  we  now  examine  them  side- 
by-side  to  assess  their  relative  strengths  and  weaknesses.  As  stated  in  the  introdnction 
(Chapter  I),  there  are  fonr  criteria  we  desire  for  an  NRBC:  speed,  accnracy,  stability, 
and  ease  of  implementation.  We  consider  each  criterion  individnally  and  compare  the 
three  NRBC  methods  against  each  other. 

Speed.  We  do  not  have  speed  comparisons  of  the  implementations,  as  a  simple 
“execntion  time”  metric  is  dependent  on  the  efficiency  of  the  code  and  the  inventive¬ 
ness  of  the  programmer.  However,  we  can  make  some  ballpark  estimates  based  on 
the  approximate  operator  connts  reqnired.  As  noted  in  their  respective  chapters,  the 
basic  Higdon  scheme  (as  antomated  by  the  Givoli-Neta  algorithm)  reqnires  0(3'^) 
operations,  bnt  the  simplihcation  of  setting  all  the  cj  to  the  same  valne  rednces  the 
operation  connt  to  O(J^).  The  Givoli-Neta  and  Hagstrom-Warbnrton  anxiliary  vari¬ 
able  methods  each  reqnire  0{J)  operations,  except  for  the  Givoli-Neta  NRBC  in  the 
presence  of  gravity,  which  reqnires  O(J^). 

Accuracy.  By  nsing  the  same  nnmerical  example  thronghont  this  disser¬ 
tation,  we  can  easily  compare  the  relative  accnracy  of  the  three  implementations. 
The  only  change  is  that  the  Higdon  scheme’s  no-advection  examples  nsed  a  semi- 
inhnite  channel  instead  of  the  infinite  channel  nsed  for  the  Givoli-Neta  and  Hagstrom- 
Warbnrton  schemes.  Even  if  we  assnme  the  error  norms  for  the  Higdon  scheme  wonld 
be  donbled  by  having  two  open  bonndaries,  we  still  see  that  the  Higdon  NRBC  error 
norms  are  approximately  80%  lower  than  Hagstrom-Warbnrton  NRBC  error  norms, 
which  in  tnrn  are  approximately  40%  lower  than  the  Givoli-Neta  NRBC  error  norms 
(for  the  basic  system,  bnt  only  10%  lower  for  the  Coriolis-  and  gravity-inflnenced 
systems).  The  only  exception  to  this  case  is  the  Higdon  scheme  with  advection  and 
Coriolis.  In  all  other  cases,  Hagstrom-Warbnrton  is  approximately  the  same  or  slightly 
better  than  Givoli-Neta,  bnt  Higdon  is  signihcanly  better  than  either  of  them.  Again, 
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this  is  a  consequence  of  the  characteristic-based  boundary  method. 

Stability.  When  we  consider  stability,  our  comparison  depends  on  the  con- 
hguration.  Comparing  the  Higdon  scheme  to  the  Givoli-Neta  scheme,  we  see  that 
Higdon  is  more  stable  when  there  is  no  advection,  but  Givoli-Neta  is  more  stable  in 
the  presence  of  non-zero  advection.  In  fact,  the  Higdon  scheme  is  completely  unstable 
in  the  non-zero  advection  case  with  Coriolis  forces;  however,  the  Givoli-Neta  scheme 
shows  some  medium-term  stability.  The  Hagstrom-Warburton  scheme  is  the  most 
stable,  although  it  and  the  Givoli-Neta  schemes  are  less  stable  in  an  open  half-plane 
with  non-zero  advection.  It  also  turns  out  that  only  the  Hagstrom-Warburton  scheme 
exhibited  any  appreciable  long-term  stability  (specihcally,  in  the  inhnite  channel  do¬ 
main),  while  the  other  two  methods  failed  for  nearly  all  conhgurations. 

Ease  of  implementation.  All  three  methods  are  fairly  straightforward  and 
easy  to  implement.  Since  all  three  methods  are  explicit,  they  can  be  computed  in 
a  separate  calculation  after  the  interior  scheme’s  calculations.  The  Higdon  scheme 
requires  a  summation  algorithm  to  automate  the  high-order  hnite-difference  approx¬ 
imations,  and  the  Givoli-Neta  and  Hagstrom-Warburton  schemes  require  a  careful 
consideration  of  characteristic  and  some  matrix  calculations.  These  drawbacks  are 
not  signihcant.  However,  Higdon  is  limited  by  the  number  of  points  in  the  domain,  and 
it  is  not  yet  possible  to  implement  Hagstrom-Warburton  on  an  open  top  boundary  in 
the  presence  of  gravity.  Furthermore,  Higdon  requires  a  careful  choice  of  interior  dis¬ 
cretization  schemes  [39,  19],  while  experiments  with  the  Hagstrom-Warburton  scheme 
for  the  scalar  wave  equation  [55]  show  that  higher-order  interior  schemes  can  be  used 
(an  option  we  have  not  attempted  to  utilize  here).  Such  high-order  schemes  may  also 
be  permissible  for  the  Givoli-Neta  method. 

Summary.  The  Higdon  scheme  has  an  enormous  advantage  in  accuracy, 
Hagstrom-Warburton  is  the  most  stable,  Givoli-Neta  and  Hagstrom-Warburton  are 
approximately  equally  fast,  and  all  three  are  equally  easy  to  implement. 
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IX.  CONCLUSIONS  AND  AREAS  FOR 
FURTHER  RESEARCH 

We  have  shown  through  the  preceding  chapters  that  high-order  non-reflecting 
boundary  conditions  can  be  applied  to  the  2-D  linearized  Euler  equations  in  a  wide 
variety  of  configurations.  The  Higdon,  Givoli-Neta,  and  Hagstrom-Warburton  tech¬ 
niques  can  all  be  used,  in  channels  and  in  open  domains,  with  and  without  advec- 
tion,  with  basic  conditions  or  with  gravity  or  Coriolis,  excepting  only  the  Hagstrom- 
Warburton  method  with  gravity.  However,  despite  the  large  amount  of  work  con¬ 
tained  herein,  there  remain  many  more  extensions  and  improvements.  These  areas 
for  further  study  include  the  following: 

1.  Application  of  the  auxiliary  variable  methods  (Givoli-Neta  and  Hagstrom- 
Warburton)  to  finite  element  models 

2.  Extension  of  the  NRBC  formulations  to  account  for  evanescent  modes  as  well 
as  the  primary  waves 

3.  Analysis  of  the  stability  of  each  scheme  for  long-time  integrations,  including 
an  exploration  of  the  conjectures  enumerated  in  Chapter  VH,  and  mitigation 
of  identified  error  sources 

4.  Application  of  the  NRBCs  to  the  full  three-dimensional  system,  including  both 
gravity  and  Coriolis  simultaneously 

5.  Incorporating  the  NRBCs  into  a  nested  environment 

6.  Finding  a  means  to  incorporate  gravity  into  the  Hagstrom-Warburton  scheme 

7.  Application  of  the  NRBCs  to  the  non-linear  system,  in  two  or  three  dimensions 

8.  Extending  this  work  to  other  linear  first-order  systems,  such  as  Maxwell’s 
equations  or  the  shallow-water  equations  (as  a  first-order  system,  not  converted 
to  the  Klein-Gordon  equation  as  in  [113]  and  elsewhere) 

In  the  introduction,  we  stated  that  one  motivation  behind  this  NRBC  develop¬ 
ment  was  to  support  the  next  generation  of  atmospheric  modeling  tools.  This  research 
has  made  significant  progress  toward  that  goal.  With  a  broad-based  finite-difference 
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implementation,  the  next  step  is  to  adapt  the  implementation  for  a  spectral  element 
system.  In  addition,  the  long-term  stability  concerns  still  need  to  be  addressed.  Al¬ 
though  we  have  not  satished  our  original  purpose,  we  have  nonetheless  developed  new 
computation  tools  for  a  broad  array  wave  propagation  models. 
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APPENDIX  A.  SIMPLIFYING  THE  EULER 

EQUATIONS 


1.  INTRODUCTION 

Eq.  (11.55)  from  Sec.  II. B  is 

dtp  +  dx{pu)  +  dy{pv)  +  d^{pw)  =  0 

dt{pu)  +  d^{pu^)  +  dy{puv)  +  d:,{puw)  +  d^p  =  fpv 
dt{pv)  +  dxipuv)  +  dy{pv^)  +  d^{pvw)  +  dyp  =  -fpu  (A.l) 
dt(pw)  +  dx{puw)  +  dy{pvw)  +  dz{pw^)  +  d^p  =  -gp 
dt{pe)  +  d^{{pe  +  p)u)  +  dy{{pe  +  p)v)  +  d^{{pe  +  p)w)  =  -gpw 


This  form  can  be  simplified,  if  we  assume  the  fluid  under  consideration  is  an 
ideal  gas  {p  =  pRT).  The  simplihed  form  requires  fewer  terms  and  is  thus  easier  and 
faster  to  calculate.  We  will  consider  the  mass,  momentum,  and  energy  equations 
separately. 

2.  MASS  EQUATION 

We  begin  the  simplihcation  process  with  the  hrst  equation  of  the  set: 

dtp  +  d^{pu)  +  dy{pv)  +  d^{pw)  =  0  (A. 2) 

This  equation  needs  no  simplihcation.  However,  we  will  expand  the  three  spatial 
derivative  terms  using  the  product  rule,  since  that  form  will  appear  in  the  subsequent 
simplihcations 

dtp  +  udxp  +  vdyp  +  wdzp  +  p  {dxU  +  dyV  +  d^w)  =  0  (A. 3) 

3.  MOMENTUM  EQUATIONS 

The  next  simplihcation  can  be  taken  with  the  next  three  equations  of  the  set: 
dt{pu)  +  d^{pu^)  +  dy{puv)  +  d^{puw)  +  d^p  =  fpv 
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dt{pv)  +  d^{puv)  +  dy{pv^)  +  d^{pvw)  +  dyp  =  -fpu 
dt{pw)  +  da;{puw)  +  dy{pvw)  +  d^{pw‘^)  +  d^p  =  -gp  (A.4) 


>  = 


>  = 


=  fpv 

fpu 


We  first  use  the  product  rule  to  separate  the  equations’  components  as  follows: 

udtp  +  udx{pu)  +  udy{pv)  +  udz{pw) 

+pdtu  +  pudxU  +  pvdyU  +  pwdzU  +  dxP 

vdtp  +  vdx{pu)  +  vdy{pv)  +  vdz{pw) 

+pdtv  +  pudxV  +  pvdyV  +  pwdzV  +  dyp 

wdtp  +  wdx{pu)  +  wdy{pv)  +  wdz{pw) 

+pdtw  +  pudxW  +  pvdyW  +  pwdzW  +  dzP 

We  then  combine  terms  to  get 

u  {dtp  +  dx{pu)  +  dy{pv)  +  dz{pw)) 

+p  {dtu  +  udxU  +  vdyU  +  wdzu)  +  dxP 

V  {dtp  +  dx{pu)  +  dy{pv)  +  dz{pw)) 


-gp 


+p  {dtv  +  udxV  +  vdyV  +  wdzv)  +  dyP 


=  fpu 

=  -fpu 


(A,5) 


w  (d,p  +  d^(pu)  +  dy(fm)  +  d,{pw)) 

> 

+p  {dtw  +  udxW  +  +  wdzw)  + 

The  terms  in  the  hrst  set  of  parentheses  matches  (A. 2),  so 
Dividing  the  remaining  terms  by  p  gives 


=  -gp  (A.6) 

we  can  eliminate  them. 


dtU  +  udxU  +  vdyU  +  wdzU  H — dxP  =  fv 

P 

dtV  +  udxV  +  vdyV  +  wdzV  ^ — dyp  =  —fu 

p 

dtW  +  udxW  +  vdyW  +  wdzW  H — dzP  =  —g 


(A.7) 


4.  ENERGY  EQUATION 

Now  the  simplihcation  process  gets  more  complicated,  using  the  hnal  equation 
of  the  set: 

dt{pe)  +  dx{{pe  +  p)u)  +  dy{{pe  +  p)v)  +  dz{{pe  +  p)w)  =  -gpw  (A. 8) 
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We  begin  with  our  state  equation  and  our  definition  of  energy 


p  =  pRT 

(A.9) 

e  =  c^T  + 

Li 

(A.  10) 

R  =  Cp  —  Cy 

(A.11) 

Combining  these  terms,  we  have 

Cv  pu^  + 

2 

(A.  12) 

Substituting  this  value  into  (A. 8)  gives 

dt  (pf  + 

+d,  ((pt  +  u  +  pu) 

+dy  ((pf  +  V  +  pv^ 

+d,  ((pf  +  e^!±£|!W)  ^  +  pui^ 

Separating  the  sums  and  using  the  product  rule  gives  the  expansion 

^dtp  +  pudtu  +  pvdtv  +  pwdtw  +  Q^p 

+  %dx{pu)  +  \dx  {pu  {u^  +  +  w^))  +  dx{pu) 

+  ^dy{pv)  +  \dy  {pv  {u^  +  +  w^))  +  5y(pn) 

+  ^^2(P'M^)  +  \dz  {puj  {u^  +  +  w^))  +  dz{pw) 

Combining  the  ^  terms  and  expanding  more  pieces  with  the  product  rule  gives 

^  {dtp  +  udxP  +  vdyp  +  wdzP  +  p  {dxU  +  dyV  +  d^w)) 

+pudtu  +  pvdtV  +  pwdtW 
+  ^2+^2+^2  (5tp  +  5^(pm)  +  5j^(pn)  +  5^(pw)) 

+(pm)  {udxU  +  vdxV  +  wdxw)  '  =  —Qpw  (A.  15) 

+{pv)  {udyU  +  vdyV  +  wdyw) 

+{pw)  {udzU  +  +  wd^w) 

+udxP  +  n5yp  +  wdzP  +  p  +  dyV  + 
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The  parenthetical  term  in  the  third  line  matches  (A. 2),  so  we  can  eliminate  it  entirely. 
In  addition,  if  we  combine  the  pu,  pv,  and  pw  terms  together,  we  get 

^  {dtp  +  udxP  +  vdyp  +  wdzP  +  p  {dxU  +  dyV  +  dzw)) 

+{pu)  {dtu  +  udxU  +  vd^v  +  wdxw) 

+{pv)  {dtv  +  udyU  +  vdyV  +  wdyw)  *  =  —gpw  (A.  16) 
+{pw)  {dtw  +  udzU  +  vdzV  +  wdzw) 

+udxP  +  vdyP  +  wdzP  +  p  {dxU  +  dyV  +  dzw) 

Based  on  the  simplihed  resnlts  in  (A.  7),  we  can  replace  the  third,  fonrth,  and  hfth 
lines  of  the  above  eqnation,  resnlting  in  the  following 

^  {dtp  +  udxP  +  vdyp  +  wdzP  +  P  {dxU  +  dyV  +  dzw)) 

+(H  {-),dxP  +  fv) 

+{pv)  (-),dyp  -  fu)  >  =  -gpw,  (A.  17) 
Hpw)  (-),dzp-g) 

+udxP  +  vdyp  +  wdzP  +  P  {dxU  +  dyV  +  dzw) 

which  we  can  expand  into 

^  {dtp  +  udxP  +  vdyp  +  wdzP  +  P  {dxU  +  dyV  +  dzw)) 

—udxP  +  fpuv  —  vdyP  —  fpuv  —  wdzP  —  gpw  =  —gpw,  (A.  18) 
+udxP  +  vdyP  +  wdzP  +  P  {dxU  +  dyV  +  dzw) 

Canceling  terms  rednces  this  eqnation  to 

-A  {dtp  +  udxP  +  vdyP  +  wdzP  +  P  {dxU  +  dyV  +  dzw))  +  p  {dxU  +  dyV  +  dzw)  =  0, 

(A.19) 

which  we  can  also  write  as 


^  {dtp  +  udxP  +  vdyp  +  wdzp)  +  (^  +  l)  T  +  dyV  +  dzw) 


From  onr  dehnition  of  R  we  have 


C|;  R  Cy  T  Cy 

R^  ^  R  ^  R 


(A.20) 

(A.21) 
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so  that  multiplying  the  above  equation  by  —  gives  the  simplihed  equation 

C-v 

dtp  +  udxP  +  vdyp  +  wdzP  +  7P  {d^u  +  dyV  +  =  0,  (A.22) 

where  7  =  ^,  and  we  have  replaced  our  energy  variable  with  the  primitive  state 
variable  for  pressure. 
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APPENDIX  B.  THE  FINITE  DIFFERENCE 
INTERIOR  SCHEME 


Let  us  briefly  analyze  the  interior  scheme  used  in  Chapters  IV~VI.  We  claim 
that  this  discretization  scheme  is  0{Sx‘^  +  Is  that  claim  true?  To  test  this 

scheme,  we  contrive  an  analytic  solution  of  sines  and  cosines 


p  =  Po  cos{kxx)  cos{kyy)  cos(u;t) 
u  =  uq  sm{kxx)  cos{kyy)  cos{ujt) 

(B.l) 

V  =  vq  cos{kxx)  sm{kyy)  cos(a;t) 
p  =  Po  cos{kxx)  cos{kyy)  cos(a;t) 

with  kx  =  ky  =  j  and  uj  =  +  ky.  On  a  40  x  40  m  domain,  these  values 

match  a  hard  wall  boundary  condition  on  all  four  sides.  We  apply  the  differential 
operators  of  (IV.  1)  to  this  equation,  and  set  the  results  as  the  right-hand-side  values  of 
(IV. 1).  These  terms  will  then  act  to  force  the  solution  to  remain  equal  to  (B.l)  within 
the  limits  of  the  discretization  error.  So  that  all  four  state  variables  are  approximately 
the  same  order  of  magnitude,  we  set  po  =  1;  7  =  1;  Po  =  1;  and  Uq  =  =  0.25. 

Since  the  leap-frog  method  requires  two  prior  time  steps,  we  use  the  analytic 
solution  to  set  the  values  for  the  hrst  two  time  steps.  We  then  use  the  discretization 
scheme  to  compute  the  next  time  step,  which  we  compare  to  the  analytic  solution  at 
that  same  time  step.  We  compute  the  average  absolute  error  at  each  interior  point 
for  each  state  variable  by 


^  (iV,  -  2)  {N,  -  2) 

where  (p  is  our  computed  state  variable,  and  p  is  the  analytic  solution  state  variable. 

For  the  hrst  test,  we  begin  with  5x  =  0.8,  and  we  halve  it  with  each  iteration. 
We  set  5y  =  0.00625  for  all  iterations,  so  that  the  discretization  error  in  y  does  not 
overwhelm  the  ^-discretization  error  we  wish  to  measure,  and  we  set  5t  =  2“®,  which 
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Sx 

Sy 

St 

Ep,  Ep 

Eu 

Ex 

0.8 

0.00625 

0.015625 

2.7983  xl0~^ 

8.0102  xl0~^ 

4.9978  xlO-^ 

0.4 

0.00625 

0.015625 

7.1  xlO"^ 

2.0223  xlO-^ 

1.2617  xlO-5 

0.2 

0.00625 

0.015625 

1.7814  xlO-5 

5.062  xlO-^ 

3.1584  xlO-® 

0.1 

0.00625 

0.015625 

4.4562  xlO-6 

1.2654  xlO-® 

7.91  xl0~^ 

Table  XL.  Discretization  errors  for  different  grid  spacings  Sx 


Sx 

Sy 

St 

Ep,  Ep 

Eu 

Eu 

0.00625 

0.8 

0.015625 

2.7983  xl0~^ 

4.9978  xl0^5 

8.0102  xlO-4 

0.00625 

0.4 

0.015625 

7.1  xlO-^ 

1.2617  xlO"® 

2.0223  xlO-^ 

0.00625 

0.2 

0.015625 

1.7814  xlO-® 

3.1584  xlO-6 

5.062  xlO-® 

0.00625 

0.1 

0.015625 

4.4562  xlO-6 

7.91  xlO-^ 

1.2654  xlO-® 

Table  XLI.  Discretization  errors  for  different  grid  spacings  Sy 

is  well  below  the  CFL  limit  for  these  conditions.  Table  XL  shows  the  discretization 
errors  for  each  state  variable  with  each  halving  of  Sx.  The  error  norms  decrease  by 
a  factor  of  almost  fonr  with  each  halving  of  Sx,  exacly  as  expected  for  a  second- 
order  method.  The  only  exception  to  this  decrease  is  the  error  norm  for  v.  We  note, 
however,  that  the  eqnation  for  v  contains  two  dy  terms  and  only  one  dx  term,  so  it 
is  not  snrprising  that  redncing  the  error  for  dx  without  simultaneously  reducing  the 
error  for  dy  wonld  have  only  a  small  impact. 

For  onr  second  test,  we  set  Sx  =  0.00625,  and  we  begin  with  Sy  =  0.8,  halving 
it  with  each  iteration.  Keeping  St  the  same  as  before,  we  get  the  resnlts  shown  in 
Table  XLI.  Note  that  these  errors  are  almost  identical  to  those  for  the  Sx  test,  with 
the  error  norms  for  u  and  v  reversed.  Hence  we  conclnde  that  the  method  is  also 
second-order  in  the  y  direction. 

As  a  third  test,  we  start  with  Sx  =  Sy  =  0.8  and  halve  both  grid  spacings  at 
each  iteration.  Using  a  larger  St  =  2~^,  we  get  the  resnlts  shown  in  Table  XLII,  and 
we  see  second-order  resnlts  in  all  fonr  state  variables. 

Finally,  we  test  St.  This  time,  we  set  Sx  =  Sy  =  0.4,  and  we  begin  with 
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Sx 

Sy 

St 

Ep,  Ep 

Eu 

E^ 

0.8 

0.8 

0.03125 

9.5561  xlO-^ 

1.5883  xlO-^ 

1.5883  xlO-^ 

0.4 

0.4 

0.03125 

2.4401  xlO-^ 

4.031  xlO-^ 

4.031  xlO-^ 

0.2 

0.2 

0.03125 

6.1292  xlO-^ 

1.012  xlO-^ 

1.012  xlO-^ 

0.1 

0.1 

0.03125 

1.5216  xlO-5 

2.5368  xlO-® 

2.5368  xlO-® 

Table  XLII.  Discretization  errors  for  different  grid  spacings  Sx  and  Sy 


Sx 

Sy 

St 

Ep,  Ep 

Eu 

Eu 

0.8 

0.8 

0.25 

0.0066075 

0.012469 

0.012469 

0.8 

0.8 

0.125 

0.0037399 

0.0063109 

0.0063109 

0.8 

0.8 

0.0625 

0.0019051 

0.0031717 

0.0031717 

0.8 

0.8 

0.03125 

0.00095561 

0.0015883 

0.0015883 

Table  XLIII.  Discretization  errors  for  different  time  steps  St 


St  =  0.25,  halving  it  with  each  iteration.  The  resnlts  are  somewhat  snrprising.  As 
shown  in  Table  XLIII,  the  error  norms  only  decrease  by  a  factor  of  two  with  each 
halving  of  St,  which  implies  a  hrst-order  method  rather  than  second-order.  The 
discretization  scheme  is  the  same.  How  can  a  scheme  which  is  second-order  in  space 
be  only  hrst-order  in  time?  Perhaps  the  difference  comes  from  how  the  discretization 
scheme  is  nsed.  For  the  spatial  derivative  approximations,  we  nse  the  known  node 
valnes  to  approximate  the  derivative,  that  is. 


dxU  ^ 


u. 


i+l 


u: 


-1 


2Sx 


(B.3) 


which  is  a  second-order  approximation,  easily  demonstrated  via  Taylor  series  expan¬ 
sion.  For  the  time  derivative,  we  approximate  it  nsing  the  eqnation  system,  and  then 
nse  that  approximation  and  the  earlier  node  valne  to  compnte  the  new  node  valne: 


2St 


dtu  =  RHS((p) 

i  dtu 


,n+l 


+2(Si9,M 


(B.4) 

(B.5) 

(B.6) 
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If  we  rewrite  the  superscript  for  the  ^  term,  we  see  that  this  method  is  in  fact 
Euler’s  method  over  25t: 

~  +  28tdtu  (B.7) 

Since  Euler’s  method  is  only  a  hrst-order  approximation,  our  time  marching  method 
is  only  hrst-order,  even  though  it  is  dehned  using  the  same  discretization  scheme 
for  our  second-order  spatial  derivative  approximations.  Thus,  it  appears  that  our 
discretization  scheme  is  in  fact  0{8x^  +  81/“^  +  8t). 

This  analysis  “feels”  wrong,  and  it  caused  some  strenuous  debates  with  the 
dissertation  supervisors.  However,  it  hts  the  observations.  Furthermore,  other  ex¬ 
periments  designed  to  test  the  leap-frog  scheme’s  performance  also  showed  similarly 
inexplicable  results.  Testing  it  as  an  ODE  solver  resulted  in  0{8x^)  performance  for 
a  single-step  and  a  global  0{8x‘^)  convergence.  Contrariwise,  a  test  using  an  analytic 
solution  to  the  heat  equation 


OtU  T  8)yyXl 

initially  showed  0{8x^f‘^)  single-step  convergence,  but  as  the  time  step  became  smaller, 
the  improvement  shrank  to  0{8x)  and  remained  there.  The  cause  of  this  performance 
is  under  investigation. 
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APPENDIX  C.  WAVE-LIKE  SOLUTIONS  OF 
THE  LINEARIZED  EULER  EQUATIONS 


We  now  derive  the  wave-like  solution  used  in  Lemma  VL2.  The  derivation  is 
similar  to  Hu’s  [71],  where  he  uses  a  non-dimensionalized  system  and  split  variables 
to  derive  the  PML  equations.  We  hrst  assume  the  existence  of  a  wave-like  solution, 
then  we  derive  the  amplitude  of  each  state  variable.  By  avoiding  contradictions,  we 
prove  that  such  a  wave-like  solution  exists.  This  wave-like  solution  has  the  form 

^^^^ikx+ily-i^t  ^  (C.l) 


where  (p  denotes  our  state  variables,  and  denotes  the  amplitudes  of  each  compo¬ 
nent.  If  we  apply  (C.l)  to  (IV. 1)  and  cancel  out  the  common  exponential  term,  we 
get 

—iojp*  +  ikuop*  +  ilvop*  +  ikpou*  +  ilpov*  =  0 
—iu)u*  +  ikuou*  +  UvqU*  +  ^p*  =  0 

(C.2) 

—iu!v*  +  ikuov*  +  Uvqv*  +  ^p*  =  0 
—iujp*  -\-  ikuop*  +  ilvop*  +  ik'ypou*  +  iljp^v*  =  0  . 

Combining  terms,  we  get 


{kuo  +  Ivq  —  ou)  p*  +  kpQU*  +  Ipov*  =  0 
{kuo  -f  Ivq  —  oj)u*  +  Vp*  =  0 
{kuo  +  Ivo  -uj)v*  +  Vp*  =  0 
k'fpou*  -1-  I'jpov*  +  {kuo  +  Ivq  —  u})p*  =  0  . 


For  acoustic  waves,  kuo  +  Ivq  —  co  0.  {kuQ  +  Ivq  —  u)  =  Q  for  entropy  and  vorticity 
waves;  see  [71]).  We  can  easily  solve  for  p*  and  p*  in  terms  of  u*  and  v*: 


Po  {kul  +  Iv*) 
oj  —  kuo  —  Ivq 
7Po  {kul  +  Iv*) 
OJ  —  kuQ  —  Ivq 


(C.4) 
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From  (C.3b,c)  we  have 


kuo  +  lvo  —  CO  ,  kuo  +  lvo  —  CO  , 
- : - u  =  - : - V 


(C.5) 


which  leads  directly  to 

u*  _  k 

V*  I 

If  we  insert  our  solution  for  p*  into  (C.3b,c),  we  get 

*  kp* 


Po 

—  kuo 

- 

■  ivo) 

k^Po {ku* 

+ 

Iv*) 

Po 

—  kuo 

- 

■  Ivof 

IPo 

k^u* 

+ 

klv* 

Po  (a;  —  kuo 

-  Ivof 

Ip* 

Po 

—  kuo 

- 

■  Ivo) 

Ijpo  {ku* 

+ 

Iv*) 

Po 

—  kuo 

- 

■  Ivof 

7Po 

klu* 

+ 

fv* 

Po  {co  —  kuo  —  Ivof 


(C.6) 


(C.7) 


(C.8) 


If  we  use  (C.6)  to  remove  I  from  the  numerator  of  (C.7a)  and  k  from  the  numerator 
of  (C.7b),  we  get 


7Po  + 

Po  {uo  —  kuQ  —  IvoY 

IPO 


u  = 


V  = 


Po  (co  —  kuQ  —  Ivof  ’ 
which,  after  a  little  algebra,  can  be  reformulated  as 


(C.9) 


_^co  —  kuo  —  Ivo  u* 

Co 

_^co  —  kuo  —  Ivq  V* 
Co  y/u*'^  + 


(C.IO) 
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Now,  let  u*  =  U  cos  9,  v*  =  U  sin  9,  for  some  U  &  C  and  9  G  [0, 2%]  measnred 
connterclockwise  from  the  positive  x  axis.  After  some  manipnlation,  we  can  explicitly 
solve  for  k  and  1: 


U!  COS  9 

Co  +  Uq  cos  9  +  Vq  sin  9 
u)  sin  9 

Co  +  Uq  cos  6*  +  no  sin  6*  ’ 


(C.ll) 


where  we  choose  the  positive  for  each  ±  in  (C.IO),  denoting  propagation  in  the  positive 
X  and  y  directions.  We  then  insert  these  valnes  of  k  and  I  into  the  eqnations  for  0* 
and  simplify,  yielding 

p*  = 

u*  =  U  cos  9 

(C.12) 

V*  =  U  sin  9 

p*  = 

^  Co 

Thns  onr  wave  solntion  becomes 


0 


if  exp 


cos  X  +  sin  y 


Co  +  Uq  cos  6*  +  no  sin  9 


(C.13) 


with  (0*  dehned  by  (C.12).  A  little  tedions  algebra  shows  that  this  solntion  satishes 
(IV. 1).  If  we  insert  this  solntion  into  (VI. 16),  we  get  the  wave  speed 


=  (co  +  Uq  cos  6*  +  no  sin  9)'^  . 


(C.14) 


Similarly,  when  we  consider  the  opposite-sign  terms  of  (C.IO),  we  again  get  a  wave-like 
solntion,  this  time  with 

cIj  =  {uq  cos  9  +  Vq  sin  9  —  cq)'^  .  (C.15) 


Since  we  do  not  know  in  advance  the  propagation  angle  of  these  plane  waves,  for 
the  pnrpose  of  dehning  the  NRBC  wave  speed,  we  assnme  the  angle  is  normal  to 
the  bonndary.  With  this  assnmption,  and  the  reqnirement  that  c^i  >  0,  we  have  the 
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following: 


—  Co  +  Mo 

Tn 

Cyj 

=  Co  +  Mo 

Tw 

Cyj 

=  Co  —  Mo 

Ts 

Cyj 

o 

1 

o 

(C.16) 
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APPENDIX  D.  DISCRETE  REFLECTION 
COEFFICIENT^A  PRELIMINARY  ANALYSIS 


Let  us  briefly  look  at  this  idea  of  quantifying  the  reflection  coefficient  for  the 
discrete  case  (see  Chapter  VII).  This  analysis  differs  from  Sec.  8.1.5  of  [24]  in  that 
we  consider  the  physical  waves,  not  merely  the  high-frequency  non-physical  compu¬ 
tational  waves  generated  by  the  hnite-difference  scheme.  For  our  initial  analysis,  we 
consider  a  Sommerfeld  (hrst-order  Higdon)  boundary  condition  in  a  1-D  domain.  A 
full  analysis  is  outside  the  scope  of  this  dissertation.  We  begin  the  work  here  to  show 
its  potential  for  future  research. 


1.  DERIVATION 

As  with  the  analysis  of  the  reflection  coefficient  for  the  continuous  equation 
(Sec.  III.B.l.a),  we  begin  by  considering  a  wave  of  the  form 

n(a;,t)  (D.l) 

traveling  left  to  right  at  unknown  speed  Cx-  (For  a  two-dimensional  wave,  we  consider 
only  the  portion  of  the  wave  traveling  parallel  to  the  x  axis.)  If  we  use  a  Sommerfeld 
condition  {dt  +  Codx)  u  =  0,  then  the  computed  solution  for  u  will  include  a  reflected 
wave  of  the  correct  magnitude  to  satisfy  the  boundary  condition.  Thus, 

u{x,  t)  =  (D.2) 


Discretizing  the  Sommerfeld  condition  with  a  backward  difference  in  x  and  t,  we  have 


u 


N 


U 


n—1 

N 


U 


+  Cq- 


N 


U 


N-l 


h 


=  0, 


(D.3) 


where  the  subscripts  N  and  N  —  1  denote  the  point  on  the  right  boundary  and  the 
point  immediately  to  its  left,  respectively;  the  superscripts  n  and  n  —  1  denote  the 
current  and  previous  time  steps,  respectively;  k  denotes  the  time  step  size;  and  h 
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denotes  the  spatial  step  size.  Thus,  at  the  right  boundary  at  time  step  n,  we  have 
X  =  Nh  and  t  =  nk.  Applying  this  discretization  to  our  wave  (D.2)  at  this  point,  we 
have 

^i{Nh—Cxnk)  _j_  h+c^nk)  _  ^i(Nh-Cx(n-l)k)  _  J^^i{Nh+Cx{n-l)k) 

+  A  _j_  j^^i{Nh+c^nk)  _  ^i({N-l)h-Cxnk)  _  J^^i{{N-l)h-c^nk)\ 


=  0,  (D.4) 


where  A  =  ^.  After  we  cancel  and  combine  terms,  we  solve  for  R  and  get 

1  -  +  A  (l  - 

^2icxnk  _  ft-icxk  ^  (^•^) 

We  are  interested  in  the  magnitude  of  R  more  than  its  actual  value.  Solving  for  this 
magnitude,  we  have 


R 


R 


1  -  +  A  (l  - 

1  —  +  A  (1  —  e“*^) 

^ _ sm{cxk) _ 

2i  (1  -  +  A  (1  -  e-*^)) 


1  - 


sm  c. 


A) 


—  icxk 

e  2  sm 


Cxk 

2 


T  Ae 


-ih 

2 


sm 


(I) 


(D.6) 


2.  IMPLICATIONS  AND  SPECULATION 

Looking  at  (D.6),  we  see  that  certain  combinations  of  c^k  and  h  could  make 
the  fractional  term  negative,  resulting  in  a  reflection  coefficient  R  with  magnitude 
greater  than  one.  In  fact,  if  we  use  the  following  constants, 

m 

Co  =  343 — 
s 

h  =  100  m 
k  =  0.18  s  , 

which  are  close  to  those  used  in  our  numerical  examples  in  this  dissertation,  we  get 
several  large  reflection  coefficients  as  Cx  varies  from  1  up  to  cq.  The  top  half  of 
Fig.  47  shows  these  reflection  coefficients.  Each  value  of  R  larger  than  one  represents 
a  potentially  unstable  wave  mode. 


158 


Constants  for  mesoscale  Euler  equations 
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Constants  for  small-scale  Klein-Gordon  equation 
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Figure  47.  Discrete  reflection  coefficients  for  varying  wave  speeds  Cx-  The  x-axis  is 
the  wave  speed  Cx'-,  the  y-axis  is  the  magnitude  of  the  reflection  coefficient  R  computed 
by  (D.6).  (Top)  Discrete  reflection  coefficients  using  constants  approximately  equal 
to  those  in  this  dissertation’s  numerical  examples  for  a  mesoscale  model.  (Bottom) 
Discrete  reflection  coefficients  using  the  same  constants  as  the  numerical  example  of 
[90]  for  the  Klein-Gordon  equation  in  a  small-scale  model. 

On  the  other  hand,  if  we  use  the  values  from  the  numerical  example  in  [90], 

^m 

s 

=  0.25  m 
=  0.125  s  , 

then  our  discrete  reflection  coefficients  are  all  less  than  one,  even  if  Cx  is  twice  as  large 
as  cq.  These  coefficients  are  plotted  in  the  bottom  half  of  Fig.  47.  These  plots  imply 
the  stability  of  the  example  in  [90]  and  the  instability  of  our  examples. 

The  problem  appears  to  be  the  scale  of  the  domain.  Perhaps  the  introduction 
of  a  scaling  factor  can  improve  the  stability.  Future  research  will  expand  this  analysis, 
to  determine  the  general  formula  for  the  Higdon  NRBC  of  order  J,  and  also  to  con¬ 
sider  the  Givoli-Neta  and  Hagstrom-Warburton  NRBGs,  after  converting  the  normal 
derivatives  to  tangential.  Ideally,  this  analysis  will  uncover  the  choice  (or  combina- 


co 

h 

k 
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tion  of  choices)  of  Cj  to  choose  to  keep  ||i?||  <  1  for  all  possible  c^.  In  the  meantime, 
this  qnick  test  has  nncovered  a  possible  reason  for  the  instabilities  exhibited  by  these 
NRBCs. 
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